C
C *** LAST REVISED ON 12-OCT-1987 17:10:02.28
C *** SOURCE FILE: [DL.GRAPHICS.LONGLIB]AUXLIB.FOR
C
C *************************************************************************
C
C	SECTION 2 OF LONGLIB GRAPHICS LIBRARY PLOT ROUTINES.
C	THESE ROUTINES ARE DEVICE INDEPENDENT ROUTINES.
C
C	THIS FILE IS COMPATIBLE WITH THE ANSI FORTRAN 77 SPECIFICATION
C	WITH THE FOLLOWING EXCEPTIONS:
C		1. TABS (^I) ARE USED TO INDENT LINES
C		2. INTEGER*2 IS USED TO SAVE SPACE IN "SYMBOL"
C		3. "PXPCGT" USES A MACHINE DEPENDENT METHOD OF
C		   EXTRACT ASCII CHARACTERS FROM AN ARRAY.
C		4. ENCODE STATEMENT AND VARIABLE FORMAT QUALIFIERS
C		   ARE USED IN "NUMBER".  ALSO, NUMBER USES A MACHINE
C		   DEPENDENT METHOD OF CHECKING IS A PASSED NUMBER IS
C		   AN INTEGER OR REAL NUMBER. THIS ROUTINE CAN BE REPLACED
C		   WITH "OLDNUMB" TO ACHIEVE NEARLY THE SAME
C		   FUNCTIONALITY.  AXIS3 ALSO USES SIMILAR APPROACH
C		   TO DECODE INPUT FORMAT.
C		5. COMMENT LINES EXTEND FROM SOURCE LINES WITH THE
C		   SEPARATOR !.
C
C *************************************************************************
C
	SUBROUTINE AXIS(X0,Y0,A0,N0,S0,T0,C0,D0,NM,ML,ICOL)
C
C	EXTENSIVELY MODIFIED BY D.LONG   7-OCT-83	AT JPL
C	REVISED 13-AUG-1987 DGL
C	+ IMPROVED APPEARANCE OF EXPONENT, ADDED EXPONENT COLOR
C
C	X0	X COORDINATE OF START OF AXIS
C	Y0	Y COORDINATE OF START OF AXIS
C	A0	CHARACTER STRING TO DESCRIBE AXIS
C	N0	NUMBER OF CHARACTERS IN STRING
C		- ON CLOCKWISE SIDE OF AXIS (NORMAL FOR X)
C		+ ON COUNTER CLOCKWISE SIDE OF AXIS (NORMAL FOR Y)
C		HUNDREDS DIGIT	= 1 NO LABELING OF AXIS--TICKS AND LINE ONLY
C		THOUSANDS DIGIT = 1 HORIZONTAL NUMBERS IN Y AXIS LABEL
C		10 THOUSANDS DIGIT=1 USE NM,ML PARAMETERS
C		100 THOUSANDS DIGIT=1 USE COLOR PARAMTERS
C	S0	LENGTH OF AXIS
C		< 0	TICKS PLACED ON OPPOSITE SIDE OF AXIS LABEL
C		= 0	NO ACTION
C		> 0	NORMAL
C	T0	ANGLE OF AXIS TO X AXIS OF PAPER
C		0.0 FOR X-AXIS
C		90.0 FOR Y-AXIS
C	C0	COORDINATE OF MINIMUM TICK ON AXIS
C	D0	SCALING DISTANCE BETWEEN "1" TICKS
C		NOTE:  THE FOLLOWING ARE REQUIRED ONLY IF ABS(NO)>10000
C	NM	NUMBER OF MINOR AXIS TICKS BETWEEN MAJOR TICKS (DEFAULT=0)
C	ML	NUMBER OF MAJOR AXIS TICKS (DEFAULT= 1 TICK/INCH)
C	ICOL	COLOR ARRAY (NEEDED ONLY IF ABS(NO)>100000
C		 ICOL(1) AXIS COLOR
C		 ICOL(2) NUMBER COLOR
C		 ICOL(3) LABEL COLOR
C		 ICOL(4) EXPONENT COLOR
C		PEN COLOR ON RETURN DEPENDS ON LAST ITEM PLOTTED
C		IN THE SEQUENCE INDICATED
C
	REAL A0(1)
	INTEGER ICOL(4)
	LOGICAL VERT,TICKS,COLOR
	CS=.15			! CHARACTER SIZE
	IF (S0.EQ.0.0) GOTO 200	! ZERO LENGTH AXIS
	VERT=.FALSE.		! NO VERTICAL NUMBERS ON HORIZONTAL AXIS
	TICKS=.TRUE.		! PUT ON TICKS
	HOR=T0
	T5=0.1			! TICK LENGTH
	B7=T5+.08		! NUMBER DISTANCE FROM AXIS
	B6=B7
	B8=0.0
	NM1=0				! NUMBER MINOR TICKS
	N2=(ABS(S0)+0.5)		! NUMBER OF MAJOR TICKS
	S1=FLOAT(N2)
	XL=1.				! INCREMENT BETWEEN MAJOR TICKS
	N1=IABS(N0)
	COLOR=.FALSE.
	IF (N1.GE.100000) THEN
		N1=MOD(N1,100000)	! USE COLOR ARRAY
		COLOR=.TRUE.
	ENDIF
	IF (N1.GE.10000) THEN
		N1=MOD(N1,10000)
		N2=IABS(ML)		! NUMBER MAJOR TICKS
		IF (N2.EQ.0) N2=1
		S1=ABS(S0)
		XL=ABS(S0)/FLOAT(N2)	! SPACING MAJOR TICKS
		NM1=IABS(NM)+1		! NUMBER MINOR TICKS
		XM=XL/FLOAT(NM1)	! SPACING MINOR TICKS
	ENDIF
	IF(N0.LT.0)GOTO 10
	B3=CS*3.8			! COUNTER-CLOCKWISE LABELING
	B4=CS+0.08
	T2=T0
	GOTO 20
10	B3=-CS*4.			! CLOCKWISE LABELING
	B4=-T5-CS-.05
	T2=T0
	T5=-T5
20	IF (N1.GE.1000)THEN
		N1=MOD(N1,1000) 	! VERTICAL NUMBERS ON HORIZONTAL AXIS
		VERT=.TRUE.
		HOR=0.0
		B4=(ABS(T5)*(1.+SIGN(1.,S0))/2.+.1)*SIGN(1.,FLOAT(N0))
		B6=.49*CS
	ENDIF
	IF (N1.GE.100) THEN
		N1=MOD(N1,100)		! NO TICKS
		TICKS=.FALSE.
	ENDIF
	T5=T5*SIGN(1.,S0)
	T1=T0*0.017453294
	T3=COS(T1)
	T4=SIN(T1)
C
	T6=T5*T3
	T5=T5*T4
	X1=X0
	Y1=Y0
	IF (COLOR) CALL PLOT(FLOAT(ICOL(1)),0.,0) ! COLOR
	DO 100 I=1,N2				! MAJOR TICKS
	IF (NM1.EQ.0) GOTO 106
		DO 105 K=1,NM1			! DO MINOR TICKS
		X2=X1+T3*FLOAT(K-1)*XM
		Y2=Y1+T4*FLOAT(K-1)*XM
		X3=X2-T5*.5
		Y3=Y2+T6*.5
		CALL PLOT(X2,Y2,3)
105		CALL PLOT(X3,Y3,2)
106	X2=X1-T5
	Y2=Y1+T6
	CALL PLOT(X2,Y2,3)
	CALL PLOT(X1,Y1,2)
	X1=X1+T3*XL
	Y1=Y1+T4*XL
	IF (T0.EQ.90.) X1=X0
100	CALL PLOT(X1,Y1,2)
	X2=X1-T5
	Y2=Y1+T6
	CALL PLOT(X2,Y2,2)			! FINISH LAST MAJOR TICK
C	CHECK FOR EXPONENT VALUE
	D1=D0					! SCALING FACTOR
	C1=C0+S1*D1				! STARTING VALUE
	E1=0.0					! EXPONENT
	IF(D1.EQ.0.0)GOTO 140
110	IF(ABS(D1).LT.10.0)GOTO 130
	D1=D1*0.1
	C1=C1*0.1
	E1=E1+1.0
	GOTO 110
120	D1=D1*10.0
	C1=C1*10.0
	E1=E1-1.0
130	IF(ABS(D1).LT.0.5)GOTO 120
140	CONTINUE			! PEN AT END OF AXIS
	IF (.NOT.TICKS) THEN
		IF (COLOR) CALL PLOT(FLOAT(ICOL(3)),0.,0) ! COLOR
		GOTO 200
	ENDIF
	IF (VERT) THEN
		C2=C1-N2*D1		! MAKE SPACE FOR VERTICAL NUMBERS
		IC=1			! ON HORIZONTAL AXIS
		IF (ABS(C2).GE.1.0) IC=IFIX(ALOG10(ABS(C2)))
		IC2=1
		IF (ABS(C1).GE.1.0) IC2=IFIX(ALOG10(ABS(C1)))
		NC1=MAX(IC,IC2)+2
		IF (C2.LT.0.0.OR.C0.LT.0.0) NC1=NC1+1
		IF (N0.GT.0.0) B4=B4+FLOAT(NC1)*CS
		B3=0.0
		B8=(.25+ABS(T5)*(SIGN(1.,S0)+1.)/2.+FLOAT(NC1)*CS)
     $			*SIGN(1.,FLOAT(N0))
	ENDIF
	X2=X1-B4*T4-B7*T3		! LOCATE CENTER NUMBER LABELS
	Y2=Y1+B4*T3-B6*T4
	N2=N2+1
	IF (COLOR) CALL PLOT(FLOAT(ICOL(2)),0.,0) ! COLOR
	DO 150 I=1,N2			! LABEL MAJOR TICKS
		CALL NUMBER(X2,Y2,CS,C1,HOR,1,-1)
		C1=C1-D1*S1/FLOAT(N2-1)
		X2=X2-T3*XL
		Y2=Y2-T4*XL
150	CONTINUE
	IF (COLOR) CALL PLOT(FLOAT(ICOL(3)),0.,0) ! COLOR
	IF (N1.GT.0) THEN		! ADD TITLE
		C2=0.0
		Y2=0.0
		CALL SYMBOL(C2,Y2,CS,A0,0.0,N1,-3) ! TITLE LENGTH
		B1=0.5*(ABS(S0)-C2)		! CENTER TITLE
		IF (E1.NE.0.) B1=B1-CS*3.	! PUT ON EXPONENT
		X2=X0+B1*T3-B3*T4-B8*T4
		Y2=Y0+B1*T4+B3*T3
		CALL SYMBOL(X2,Y2,CS,A0,T2,N1,-1)
	ELSE
		C2=0.0
		B1=0.5*ABS(S0)
		X2=X0+B1*T3-B3*T4-B8*T4
		Y2=Y0+B1*T4+B3*T3
	ENDIF
	IF (E1.EQ.0.0) GOTO 200		! NO EXPONENT
	IF (COLOR) CALL PLOT(FLOAT(ICOL(4)),0.,0) ! COLOR
	C2=C2+CS
	X2=X2+C2*T3
	Y2=Y2+C2*T4
	CALL SYMBOL(X2,Y2,CS,'(X10',T2,4,-1)
	X2=X2+CS*3.75*T3-CS*T4*0.4
	Y2=Y2+CS*3.75*T4+CS*T3*0.4
	CALL NUMBER(X2,Y2,CS,E1,T2,-1,-1)
	B2=0.8+AINT(ALOG10(ABS(E1)))
	IF (E1.LT.0.0) B2=B2+1
	X2=X2+B2*CS*T3+CS*T4*0.4
	Y2=Y2+B2*CS*T4-CS*T3*0.4
	CALL SYMBOL(X2,Y2,CS,')',T2,1,-1)
200	RETURN
	END
C
C
	SUBROUTINE AXIS2(X0,Y0,A0,N0,S0,T0,C0,D0,NM,
     $		NN,ML,TS,ND,SM,ICOL)
C
C	WRITTEN BY D.LONG   17-OCT-1983		AT JPL
C	REVISED 13-AUG-1987 DGL
C	+ IMPROVED APPEARANCE OF EXPONENT, ADDED EXPONENT COLOR
C
C	THIS VERSION OF AXIS INCLUDES MINOR TICKS ON AXIS AND A MORE
C	FLEXIBLE METHOD OF SPECIFYING AXIS PARAMETERS.
C
C	X0	X COORDINATE OF START OF AXIS
C	Y0	Y COORDINATE OF START OF AXIS
C	A0	CHARACTER STRING TO DESCRIBE AXIS
C	N0	NUMBER OF CHARACTERS IN STRING
C		- ON CLOCKWISE SIDE OF AXIS (NORMAL FOR X)
C		+ ON COUNTER CLOCKWISE SIDE OF AXIS (NORMAL FOR Y)
C		HUNDREDS DIGIT	= 1 NO LABELING OF AXIS--TICKS AND LINE ONLY
C		THOUSANDS DIGIT = 1 HORIZONTAL NUMBERS IN Y AXIS LABEL
C		10 THOUSANDS DIGIT=1 USE NM,ML,ETC. PARAMETERS
C		100 THOURSANDS DIGIT=1 USE COLOR PARAMTERS
C	S0	LENGTH OF AXIS
C		< 0	TICKS PLACED ON OPPOSITE SIDE OF AXIS LABEL
C		= 0	NO ACTION
C		> 0	NORMAL
C	T0	ANGLE OF AXIS TO X AXIS OF PAPER
C		0.0 FOR X-AXIS
C		90.0 FOR Y-AXIS
C	C0	COORDINATE OF MINIMUM TICK ON AXIS
C	D0	SCALING DISTANCE BETWEEN  1" TICKS
C		NOTE:  THE FOLLOWING ARE REQUIRED ONLY IF ABS(NO)>10000
C	NM	NUMBER OF MINOR AXIS TICKS BETWEEN MAJOR TICKS (DEFAULT=0)
C	NN	NNth MINOR TICK FROM MAJOR IS HIGHLIGHTED
C		(DEFAULT=0 INDICATES NO MINOR TICKS HIGHLIGHTED)
C	ML	NUMBER OF MAJOR AXIS TICKS (DEFAULT= 1 TICK/INCH)
C		< 0 THEN USE FOLLOWING VARABLES
C	TS	CHARACTER SIZE OF TITLE AND NUMBER (DEFAULT=.15)
C		< 0 THEN DO NOT AUTO SCALE BY (x10 TO POWER)
C	ND	NUMBER OF DIGITS TO RIGHT OF DECIMAL POINTS (DEFAULT=1)
C	SM	MAJOR TICK LENGTH (DEFAULT= .1) (MINOR TICK=1/2 MAJOR TICK)
C	ICOL	COLOR ARRAY (NEEDED ONLY IF ABS(NO)>100000
C		 ICOL(1) AXIS COLOR
C		 ICOL(2) NUMBER COLOR
C		 ICOL(3) LABEL COLOR
C		 ICOL(4) EXPONENT COLOR
C		PEN COLOR ON RETURN DEPENDS ON LAST ITEM PLOTTED
C		IN THE SEQUENCE INDICATED
C
	REAL A0(1)
	INTEGER ICOL(4)
	LOGICAL VERT,TICKS,COLOR,SCALE
	CS=.15			! CHARACTER SIZE
	IF (S0.EQ.0.0) GOTO 200	! ZERO LENGTH AXIS
	VERT=.FALSE.		! NO VERTICAL NUMBERS ON HORIZONTAL AXIS
	TICKS=.TRUE.		! PUT ON TICKS
	SCALE=.TRUE.		! (x10 TO POWER SCALING)
	HOR=T0
	NDD=1			! NUMBER OF DIGITS TO RIGHT OF DECIMAL
	T5=0.1			! TICK LENGTH
	B7=T5+.08		! NUMBER DISTANCE FROM AXIS
	B6=B7
	B8=0.0
	NM1=0				! NUMBER MINOR TICKS
	N2=(ABS(S0)+0.5)		! NUMBER OF MAJOR TICKS
	S1=N2
	XL=1.				! INCREMENT BETWEEN MAJOR TICKS
	N1=IABS(N0)
	COLOR=.FALSE.
	IF (N1.GT.100000) THEN
		N1=MOD(N1,100000)	! USE COLOR ARRAY
		COLOR=.TRUE.
	ENDIF
	IF (N1.GT.10000) THEN
		N1=MOD(N1,10000)
		N2=IABS(ML)		! NUMBER MAJOR TICKS
		S1=ABS(S0)
		IF (N2.EQ.0) N2=1
		XL=ABS(S0)/FLOAT(N2)	! SPACING MAJOR TICKS
		NM1=IABS(NM)+1		! NUMBER MINOR TICKS
		XM=XL/FLOAT(NM1)	! SPACING MINOR TICKS
		IF (ML.LT.0) THEN
			CS=ABS(TS)	! DIFFERENT TITLE SIZE
			IF (CS.EQ.0.) CS=.15
			NDD=IABS(ND)
			IF (TS.LT.0) SCALE=.FALSE.	! DO NOT SCALE
			T5=ABS(SM)			! NEW TICK LENGTH
			IF (T5.EQ.0.) T5=.1
		ENDIF
	ENDIF
	IF (N1.GT.1000)THEN
		N1=MOD(N1,1000) 	! VERTICAL NUMBERS ON HORIZONTAL AXIS
		VERT=.TRUE.
		HOR=0.0
		B4=(ABS(T5)*(1.+SIGN(1.,S0))/2.+.1)*SIGN(1.,FLOAT(N0))
		B6=.49*CS
	ENDIF
	IF (N1.GT.100) THEN
		N1=MOD(N1,100)		! NO TICKS
		TICKS=.FALSE.
	ENDIF
	IF(N0.LT.0)GOTO 10
	B3=CS*(2.8+NDD) 		! COUNTER-CLOCKWISE LABELING
	B4=CS+T5
	T2=T0
	GOTO 20
10	B3=-CS*(3.+NDD) 		! CLOCKWISE LABELING
	B4=-T5-CS
	T2=T0
	T5=-T5
20	CONTINUE
	T5=T5*SIGN(1.,S0)
	T1=T0*0.017453294
	T3=COS(T1)
	T4=SIN(T1)
C
	T6=T5*T3
	T5=T5*T4
	X1=X0
	Y1=Y0
	IF (COLOR) CALL PLOT(FLOAT(ICOL(1)),0.,0) ! COLOR
	DO 100 I=1,N2				! MAJOR TICKS
	IF (NM1.EQ.0) GOTO 106
		DO 105 K=1,NM1			! DO MINOR TICKS
		X2=X1+T3*FLOAT(K-1)*XM
		Y2=Y1+T4*FLOAT(K-1)*XM
		IF (K-1.EQ.NN.AND.NN.NE.0) THEN
			HMT=.8
		ELSE
			HMT=.5
		ENDIF
		X3=X2-T5*HMT
		Y3=Y2+T6*HMT
		CALL PLOT(X2,Y2,3)
105		CALL PLOT(X3,Y3,2)
106	X2=X1-T5
	Y2=Y1+T6
	CALL PLOT(X2,Y2,3)
	CALL PLOT(X1,Y1,2)
	X1=X1+T3*XL
	Y1=Y1+T4*XL
	IF (T0.EQ.90.) X1=X0
100	CALL PLOT(X1,Y1,2)
	X2=X1-T5
	Y2=Y1+T6
	CALL PLOT(X2,Y2,2)			! FINISH LAST MAJOR TICK
C	CHECK FOR EXPONENT VALUE
	D1=D0					! SCALING FACTOR
	C1=C0+D1*S1				! STARTING VALUE
	E1=0.0					! EXPONENT
	IF (.NOT.SCALE) GOTO 140
	IF(D1.EQ.0.0) GOTO 140
110	IF(ABS(D1).LT.10.0)GOTO 130
	D1=D1*0.1
	C1=C1*0.1
	E1=E1+1.0
	GOTO 110
120	D1=D1*10.0
	C1=C1*10.0
	E1=E1-1.0
130	IF(ABS(D1).LT.0.5)GOTO 120
140	CONTINUE			! PEN AT END OF AXIS
	IF (.NOT.TICKS) THEN
		IF (COLOR) CALL PLOT(FLOAT(ICOL(3)),0.,0) ! COLOR
		GOTO 200
	ENDIF
	IF (VERT) THEN
		C2=C1-N2*D1		! MAKE SPACE FOR VERTICAL NUMBERS
		IC=1			! ON HORIZONTAL AXIS
		IF (ABS(C2).GE.1.0) IC=IFIX(ALOG10(ABS(C2)))
		IC2=1
		IF (ABS(C1).GE.1.0) IC2=IFIX(ALOG10(ABS(C1)))
		NC1=MAX(IC,IC2)+2
		IF (C2.LT.0.0.OR.C0.LT.0.0) NC1=NC1+1
		IF (N0.GT.0.0) B4=B4+FLOAT(NC1+NDD)*CS
		B3=0.0
		B8=(.25+ABS(T5)*(SIGN(1.,S0)+1.)/2.+FLOAT(NC1+NDD)*CS)
     $			*SIGN(1.,FLOAT(N0))
	ENDIF
	X2=X1-B4*T4-B7*T3		! LOCATE CENTER NUMBER LABELS
	Y2=Y1+B4*T3-B6*T4
	N2=N2+1
	IF (COLOR) CALL PLOT(FLOAT(ICOL(2)),0.,0) ! COLOR
	NDDD=NDD
	IF (NDD.EQ.0) NDDD=-1
	DO 150 I=1,N2			! LABEL MAJOR TICKS
		CALL NUMBER(X2,Y2,CS,C1,HOR,NDDD,-1)
		C1=C1-D1*S1/FLOAT(N2-1)
		X2=X2-T3*XL
		Y2=Y2-T4*XL
150	CONTINUE
	IF (N1.NE.0) THEN
		C2=0.0
		Y2=0.0
		CALL SYMBOL(C2,Y2,CS,A0,0.,N1,-3)
		B1=0.5*(ABS(S0)-C2)		! CENTER TITLE
		IF (E1.NE.0.0) B1=B1-CS*3.	! PUT ON EXPONENT SPACE
		X2=X0+B1*T3-B3*T4-B8*T4
		Y2=Y0+B1*T4+B3*T3
		IF (COLOR) CALL PLOT(FLOAT(ICOL(3)),0.,0) ! COLOR
		CALL SYMBOL(X2,Y2,CS,A0,T2,N1,-1)
	ELSE
		C2=0.0
		B1=0.5*ABS(S0)
		X2=X0+B1*T3-B3*T4-B8*T4
		Y2=Y0+B1*T4+B3*T3
	ENDIF
	IF (E1.EQ.0.0) GOTO 200		! NO EXPONENT
	IF (COLOR) CALL PLOT(FLOAT(ICOL(4)),0.,0) ! COLOR
	C2=C2+CS
	X2=X2+C2*T3
	Y2=Y2+C2*T4
	CALL SYMBOL(X2,Y2,CS,'(X10',T2,4,-1)
	X2=X2+3.75*CS*T3-CS*T4*0.4
	Y2=Y2+3.75*CS*T4+CS*T3*0.4
	CALL NUMBER(X2,Y2,CS,E1,T2,-1,-1)
	B2=0.8+AINT(ALOG10(ABS(E1)))
	IF (E1.LT.0.0) B2=B2+1
	X2=X2+B2*CS*T3+CS*T4*0.4
	Y2=Y2+B2*CS*T4-CS*T3*0.4
	CALL SYMBOL(X2,Y2,CS,')',T2,1,-1)
200	RETURN
	END
C
C
	SUBROUTINE AXIS3(X0,Y0,T,N0,S0,A0,B0,C0,D0,E0,F0,ICOL)
C
C	WRITTEN BY DGL  17-OCT-1983
C	REVISED BY DGL  28-AUG-1987
C
C	PLOTS A SINGLE COORDINATE AXIS USING SYMBOL AND NUMBER.
C	THIS VERSION OF AXIS CAN HANDLE NON-INTEGER AXIS LENGTH
C	AND PERMITS FORMATTING OF THE NUMERIC LABELS.  ROUTINE DECODES
C	THE LABEL FORMAT USING MACHINE DEPENDENT TEST FOR FORMAT
C	SPEC PASSED AS INTEGER OR FLOATING POINT -- SEE NUMBER.
C
C	X0	X COORDINATE OF START OF AXIS
C	Y0	Y COORDINATE OF START OF AXIS
C	T	CHARACTER STRING TO DESCRIBE AXIS
C	N0	NUMBER OF CHARACTERS IN STRING
C		- ON CLOCKWISE SIDE OF AXIS (NORMAL FOR X)
C		+ ON COUNTER CLOCKWISE SIDE OF AXIS (NORMAL FOR Y)
C		HUNDREDS DIGIT	= 1 NO LABELING OF AXIS--TICKS AND LINE ONLY
C		THOUSANDS DIGIT = 1 HORIZONTAL NUMBERS IN Y AXIS LABEL
C	S0	LENGTH OF AXIS IN INCHES
C		< 0	TICKS PLACED ON OPPOSITE SIDE OF AXIS LABEL
C		= 0	NO ACTION (NOP)
C		> 0	NORMAL
C	A0	ANGLE OF AXIS WITH RESPECT TO X AXIS OF PAPER
C			0.0 FOR X-AXIS
C			90.0 FOR Y-AXIS
C	B0	MINIMUM VALUE ON TICK AXIS
C	C0	MAXIMUM VALUE ON TICK AXIS
C	D0	INT(D0) = NUMBER OF MAJOR AXIS TICKS
C		INT((INT(D0)-D0)*100) = NUMBER OF MINOR AXIS TICKS BETWEEN MAJOR TICKS
C		INT(MOD((INT(D0)-D0)*10000,100) = NUMBER OF SUB MINOR AXIS TICKS 
C			BETWEEN MINOR TICKS
C	E0	CHARACTER SIZE OF TITLE AND NUMBER (IF E0=0 DEFAULTS TO .15)
C			< 0 THEN DO NOT AUTO SCALE BY (x10 TO POWER)
C	F0	NUMBER SPECIFICATION (FORMAT FX.Y)
C		 (SEE NUMBER)
C		 WHEN AN INTEGER IS PASSED AS F0, ROUTINE DETECTS
C		 THIS AND DECODES AS PER NUMBER.  IF AUTO SCALING IS
C		 ENABLED AN INPUT INTEGER OF -1 PRODUCES F0=1003.0
C		 WHILE AN INTEGER 0<=N<=12 PRODUCES F0=2+N*1.01.  WHEN
C		 AUTO SCALING IS DISABLED, FORMAT FOR NUMBER IS SET
C		 ACCORDING TO SIZE OF ACTUAL AXIS LABELING NUMBERS.
C	ICOL	OPTIONAL COLOR ARRAY (USED ONLY IF 100000'S DIGIT <> 0)
C		 ICOL(1) = AXIS LINE AND TICK COLOR
C		 ICOL(2) = NUMERIC LABEL COLOR
C		 ICOL(3) = TITLE COLOR
C		 ICOL(4) = EXPONENT COLOR
C		PEN COLOR ON RETURN WILL DEPEND ON OPTIONS IN THE
C		ORDER LISTED ABOVE
C
	INTEGER T(1),ICOL(1)
	DATA SPACE/0.08/		! MIN SPACING BETWEEN ITEMS
	EQUIVALENCE (NN,FA)		! FOR INTEGER FORMAT DECODE
	ROTX(X,Y)=CO*X-SI*Y+X01		! ROTATION MATRIX
	ROTY(X,Y)=SI*X+CO*Y+Y01
C
	TL=0.1				! TICK LENGTH
	X01=X0
	Y01=Y0
	ANG=A0				! ROTATION ANGLE
	E1=E0				! CHARACTER SIZE
	CS=ABS(E1)
	IF (CS.EQ.0.) CS=.15
	IF (S0.EQ.0.0) GOTO 1000	! ZERO LENGTH AXIS
	X1=D0*1.000002
	NJT=ABS(X1)			! NUMBER OF MAJOR TICKS
	NNT=100.0*(ABS(X1)-NJT)		! NUMBER OF MINOR TICKS
	NST=100.0*((ABS(X1)-NJT)*100.0-NNT)	! NUMBER OF SUB-MINOR TICKS
	IF (NJT.LT.2) NJT=2
	XJ=ABS(S0)/(NJT-1)		! INCREMENT BETWEEN MAJOR TICKS
	XN=XJ
	IF (NNT.NE.0) XN=XN/(NNT+1)	! INCREMENT BETWEEN MINOR TICKS
	XS=XN/(NST+1)			! INCREMENT BETWEEN SUB-MINOR TICKS
	N1=MOD(N0,100000)
	IF (IABS(N1)/1000.NE.0) HOR=90.0 ! ROTATION ANGLE
	CO=COS(ANG*0.017453294)		! AXIS ANGLE ROTATION
	SI=SIN(ANG*0.017453294)
	HOR=ANG				! ANGLE OF NUMBER LABELS
	NC=MOD(IABS(N1),100)		! NUMBER OF CHARACTERS IN TITLE
C
C	DECODE NUMBER FORMAT
C	CHECK TO SEE IF INTEGER WAS PASSED FOR F0 (MACHINE DEPENDENT)
C
	FA=F0
        IF (ABS(FA).GT.1.E30) THEN		! INPUT WAS INTEGER -1
		FA=1003.0		! DEFAULT AUTO SCALING FORMAT
		IF (E1.LT.0.0) THEN	! NO AUTO SCALING SO MAKE
		  NG=1			! FORMAT TO FIT
		  IF (B0.NE.0.0) NG=MAX(NG,INT(ALOG10(
     $					ABS(B0))+0.001))
		  IF (C0.NE.0.0) NG=MAX(NG,INT(ALOG10(
     $					ABS(C0))+0.001))
		  IF (B0.LT.0.0.OR.C0.LT.0) NG=NG+1
		  FA=1000.0+FLOAT(NG)
		ENDIF
	ENDIF
	IF (FA.EQ.0.0) THEN		! INPUT WAS INTEGER
		ND=NN			! INPUT INTEGER VALUE
		FA=3.+FLOAT(ND)*(1.01)	! DEFAULT FORMAT FOR NOT AUTOSCALE
		IF (E1.LT.0.0) THEN	! NO AUTO SCALING TO MAKE FORMAT
		  NG=2			! WHICH FITS
		  IF (B0.NE.0.0) NG=MAX(NG,INT(ALOG10(
     $					ABS(B0))+0.001)+1)
		  IF (C0.NE.0.0) NG=MAX(NG,INT(ALOG10(
     $					ABS(C0))+0.001)+1)
		  IF (B0.LT.0.0.OR.C0.LT.0) NG=NG+1
		  FA=FLOAT(NG)+FLOAT(ND)*1.01
		ENDIF
	ENDIF
	ND=MOD(ABS(FA),1000.0)		! NUMBER OF DIGITS
	NG=100.*(MOD(ABS(FA),1000.)-ND+0.0001)! NUMBER OF DIGITS RIGHT OF D.P.
	IF (NG.GT.17) NG=NG/10		! CORRECT INPUT SIMPLE ERRORS
	IF (FA.LT.0.0) ND=ND-4		! EXPONENTIAL NOTATION
	IF (ND.LE.0) ND=NG
	NDD=ND
	IF (FA.LT.0.0) NDD=ND+4		! EXPONENTIAL NOTATION
	IF (ABS(FA).GT.1000) NG=-1	! FORMATTED INTEGER
C
	TL1=TL
	IF (S0.LT.0.0) TL1=-TL		! REVERSE SIDE OF TICKS
	IF (S0.LT.0.0) TL=0.0		! REVERSE SIDE OF TICKS
	IF (IABS(N1)/1000.NE.0) GOTO 10
		DNX=-NDD*CS/2.0		! NUMBER LABEL DISTANCE FROM AXIS
		DNY=-TL-SPACE-CS	! NUMBER LABEL DISTANCE FROM AXIS
		DTY=DNY-CS-SPACE	! TITLE DISTANCE FROM AXIS
		GOTO 20
10	CONTINUE			! HORIZONTAL NUMBERS ON VERTICAL AXIS
		DNX=-CS/2.0		! NUMBER LABEL DISTANCE FROM AXIS
		DNY=-TL-SPACE		! NUMBER LABEL DISTANCE FROM AXIS
		DTY=DNY-NDD*CS-2.*SPACE	! TITLE DISTANCE FROM AXIS
		HOR=ANG-90.0
20	CONTINUE
	IF (N1.LT.0) GOTO 30		! CLOCKWISE TITLES
		DNY=-DNY-CS		! COUNTER-CLOCKWISE TITLES
		DTY=-DTY-CS
		TL1=-TL1		! CHANGE SIDES OF TICKS
		IF (IABS(N1).LT.1000) GOTO 30
			DNY=DNY+CS*NDD
			DTY=DNY+SPACE
30	X1=0.0				! FIRST MAJOR TICK
	Y1=0.0
	Y2=-TL1
	IF (IABS(N0).GE.100000) CALL PLOT(FLOAT(ICOL(1)),0.,0)
	CALL PLOT(ROTX(X1,Y1),ROTY(X1,Y1),3)
	CALL PLOT(ROTX(X1,Y2),ROTY(X1,Y2),2)
	DO 100 I=1,NJT-1		! MAJOR TICKS
		CALL PLOT(ROTX(X1,Y1),ROTY(X1,Y1),3)
		X1=XJ*I
		Y2=-TL1
		CALL PLOT(ROTX(X1,Y1),ROTY(X1,Y1),2)
		CALL PLOT(ROTX(X1,Y2),ROTY(X1,Y2),2)
		DO 110 J=1,NNT+1	! MINOR TICKS
			Y2=-TL1*0.7
			X2=X1+XN*J-XJ
			CALL PLOT(ROTX(X2,Y1),ROTY(X2,Y1),3)
			CALL PLOT(ROTX(X2,Y2),ROTY(X2,Y2),2)
			Y2=-TL1*0.4
			DO 120 K=1,NST	! SUB MINOR TICKS
				X3=X2+XS*K-XN
				CALL PLOT(ROTX(X3,Y1),ROTY(X3,Y1),3)
				CALL PLOT(ROTX(X3,Y2),ROTY(X3,Y2),2)
120			CONTINUE
110		CONTINUE
100	CONTINUE
	IF (MOD(IABS(N1),1000).GT.100) GOTO 1000 ! NO LABELING
	XS=0.0					! EXPONENT
	IF (E1.LT.0.) GOTO 140			! NO AUTO SCALING
C
C	COMPUTE AUTO EXPONENT SCALING SO THAT THE FORMATED LABEL
C	HAS THE INTEGER PORTION FILLED AS MUCH AS POSSIBLE
C
	I=ND-NG-1
	IF (B0.LT.0.0.OR.C0.LT.0.0) I=I-1
	IF (B0.NE.0.0) THEN
		X1=ALOG10(ABS(B0)+1.E-30)
		IF (X1.LT.0.0.AND.ABS(AINT(x1-0.001)-X1).
     $			GT.0.001) X1=X1-1.0
		IF (X1.GE.0.0) X1=X1+1.0
		X1=AINT(X1)
	ELSE
		X1=0.0
	ENDIF
	IF (C0.NE.0.0) THEN
		Y1=ALOG10(ABS(C0)+1.E-30)
		IF (Y1.LT.0.0.AND.ABS(AINT(Y1-0.001)-Y1).
     $			GT.0.001) Y1=Y1-1.0
		IF (Y1.GE.0.0) Y1=Y1+1.0
		Y1=AINT(Y1)
	ELSE
		Y1=0.0
	ENDIF
	X2=MIN(X1,Y1)
	X3=MAX(X1,Y1)
	IF (X3.LT.0.0) X3=X3+1
	IF (X2.LT.0.0.AND.NG.LE.2-X2) XS=ND-NG-1-X2
	IF (I.LT.X3+XS) XS=I-X3
140	Y1=DNY
	Y2=(C0-B0)/(NJT-1)
	IF (IABS(N0).GE.100000) CALL PLOT(FLOAT(ICOL(2)),0.,0)
	E1=XS					! EXPONENT VALUE
	DO 150 I=1,NJT				! LABEL MAJOR TICKS
		X1=(I-1)*XJ+DNX
		C1=(Y2*(I-1)+B0)*10.**E1
		CALL NUMBER(ROTX(X1,Y1),ROTY(X1,Y1),CS,C1,HOR,F0,-1)
150	CONTINUE
C
C	PLOT TITLE
C
	IF (NC.NE.0) THEN
	   Y1=0.
	   X1=0.
	   CALL SYMBOL(X1,Y1,CS,T,0.,NC,-3)	! GET TITLE LENGTH
	   DTX=(ABS(S0)-X1)/2.			! CENTER TITLE
	   IF (E1.NE.0.0) DTX=DTX-CS*3.0	! ADD EXPONENT SPACE
	   IF(IABS(N0).GE.100000)CALL PLOT(FLOAT(ICOL(3)),0.,0)
	   CALL SYMBOL(ROTX(DTX,DTY),ROTY(DTX,DTY),CS,T,ANG,NC,-1)
	   DTX=DTX+X1
	ELSE
	   DTX=ABS(S0)/2.
	   IF (E1.NE.0.0) DTX=DTX-CS*3.0	! ADD EXPONENT SPACE
	ENDIF
	X1=DTX+CS/2.
	Y1=DTY
	IF (E1.EQ.0.0) GOTO 1000		! NO EXPONENT
	IF(IABS(N0).GE.100000)CALL PLOT(FLOAT(ICOL(4)),0.,0)
	E1=-E1
	CALL SYMBOL(ROTX(X1,Y1),ROTY(X1,Y1),CS,'(X10',ANG,4,-1)
	X1=X1+3.75*CS
	Y1=Y1+CS*0.4
	CALL NUMBER(ROTX(X1,Y1),ROTY(X1,Y1),CS,E1,ANG,-1,-1)
	X2=AINT(ALOG10(ABS(E1)))+0.75
	IF (E1.LT.0.0) X2=X2+1.0
	X1=X1+X2*CS
	Y1=Y1-CS*0.4
	CALL SYMBOL(ROTX(X1,Y1),ROTY(X1,Y1),CS,')',ANG,1,-1)
1000	RETURN
	END
C
C
	SUBROUTINE GRID(X,Y,DX,DY,NX,NY)
C
C	ROUTINE TO PRODUCE CARTESIAN GRID
C
C	CREATED BY D. LONG    AUG, 1983	AT JPL
C
C	X,Y	CORRDINATES OF BOTTOM LEFT CORNER
C	DX,DY	SPACING OF GRID LINES IN X AND Y DIRECTIONS
C	NX,NY	NUMBER OF GRIDS IN X AND Y DIRECTIONS
C		IF NY > 0 AND NX > 0 THEN SOLID GRID DRAWN
C		IF NY < 0 OR NX < 0  TICK MARKS USED
C		IF NY < 0 AND NX < 0 THEN GRID AREA IS BOXED
C
	IF (NX.LT.0.OR.NY.LT.0) GOTO 200
	IF (NX.EQ.0.OR.NY.EQ.0) RETURN
	CALL PLOT(X,Y,3)			! SOLID GRID
	DO 100 IX=0,NX
		X1=X+IX*DX
		X2=X+(IX+1)*DX
		Y1=Y+NY*DY
		CALL PLOT(X1,Y1,2)
		CALL PLOT(X2,Y,3)
 100	CONTINUE
	CALL PLOT(X,Y,3)
	DO 150 IY=0,NY
		Y1=Y+IY*DY
		Y2=Y+(IY+1)*DY
		X1=X+NX*DX
		CALL PLOT(X1,Y1,2)
		CALL PLOT(X,Y2,3)
 150	CONTINUE
	CALL PLOT(X,Y,3)
	RETURN
 200	NX1=IABS(NX)				! TICKED GRID
	NY1=IABS(NY)
	EX=X+FLOAT(NX1)*DX
	EY=Y+FLOAT(NY1)*DY
	IF (ISIGN(1,NX).EQ.ISIGN(1,NY)) THEN	! BOX GRID AREA
		CALL PLOT(X,Y,3)
		CALL PLOT(EX,Y,2)
		CALL PLOT(EX,EY,2)
		CALL PLOT(X,EY,2)
		CALL PLOT(X,Y,2)
	ENDIF
	T=.05					! TICK SIZE
	DO 300 IX=0,NX1
		DO 300 IY=0,NY1
			X0=X+FLOAT(IX)*DX
			Y0=Y+FLOAT(IY)*DY
			X1=X0+T
			X2=X0-T
			Y1=Y0+T
			Y2=Y0-T
			IF (X2.LT.X)  X2=X
			IF (X1.GT.EX) X1=EX
			IF (Y2.LT.Y)  Y2=Y
			IF (Y1.GT.EY) Y1=EY
			CALL PLOT(X1,Y0,3)
			CALL PLOT(X2,Y0,2)
			CALL PLOT(X0,Y1,3)
			CALL PLOT(X0,Y2,2)
300	CONTINUE
	RETURN
	END
C
C
	SUBROUTINE LGRID(X,Y,DX,DY,NX,NY,IF)
C
C	ROUTINE TO PRODUCE LOGRITHMIC GRID
C
C	CREATED BY D. LONG    NOV, 1983		AT JPL
C
C	X,Y	COORDINATES OF BOTTOM LEFT CORNER
C	DX,DY	SPACING OF GRID LINES IN X AND Y DIRECTIONS
C	NX,NY	NUMBER OF GRIDS IN X AND Y DIRECTIONS
C		IF NY > 0 AND NX > 0 THEN LINEAR GRID
C		IF NY < 0 AND NX < 0 THEN LOG X AND LINEAR Y
C		IF NY < 0 AND NX < 0 THEN LOG X AND LOG Y
C		IF NY > 0 AND NX > 0 THEN LINEAR Y AND LOG Y
C	IF	OPTION FLAG
C		= 0 SOLID.  MINOR LINES FOR LOG SPACING
C		= 1 DOTTED. MINOR LINES FOR LOG SPACING
C		= 2 SOLID.  MINOR TICKS FOR LOG SPACING
C
	IF (NX.EQ.0.OR.NY.EQ.0) RETURN
	IF (IF.EQ.1) CALL NEWPEN(1)		! DOTTED LINES
	XM=ABS(NX)*DX+X
	YM=ABS(NY)*DY+Y
	INX1=IABS(NX)+1
	INY1=IABS(NY)+1
	CALL PLOT(X,Y,3)
	CALL PLOT(X,YM,2)
	DO 100 IX=1,IABS(NX)
		IF (NX.LT.0) THEN
			X0=(IX-1)*DX+X
			DO 50 IXM=2,9
			    X1=ALOG10(FLOAT(IXM))*DX+X0
			    IF (IF.GE.2) THEN
				 DO 30 IY=1,INY1
				    Y0=(IY-1)*DY+Y-.08
				    Y1=Y0+.16
				    IF (IY.EQ.1.AND.DY.GT.0.0) Y0=Y
				    IF (IY.EQ.1.AND.DY.LT.0.0) Y1=Y
				    IF (IY.EQ.INY1.AND.DY.GT.0.0) Y1=YM
				    IF (IY.EQ.INY1.AND.DY.LT.0.0) Y0=YM
				    CALL PLOT(X1,Y0,3)
				    CALL PLOT(X1,Y1,2)
30				CONTINUE
			    ELSE
				CALL PLOT(X1,Y,3)
				CALL PLOT(X1,YM,2)
			    ENDIF
50			CONTINUE
		ENDIF
	X0=IX*DX+X
	CALL PLOT(X0,Y,3)
	CALL PLOT(X0,YM,2)
100	CONTINUE
	CALL PLOT(X,Y,3)
	CALL PLOT(XM,Y,2)
	DO 200 IY=1,IABS(NY)
		IF (NY.LT.0) THEN
			Y0=(IY-1)*DY+Y
			DO 150 IYM=2,9
			    Y1=ALOG10(FLOAT(IYM))*DY+Y0
			    IF (IF.GE.2) THEN
				 DO 130 IX=1,INX1
				    X0=(IX-1)*DX+X-.08
				    X1=X0+.16
				    IF (IX.EQ.1.AND.DX.GT.0.0) X0=X
				    IF (IX.EQ.1.AND.DX.LT.0.0) X1=X
				    IF (IX.EQ.INX1.AND.DX.GT.0.0) X1=XM
				    IF (IX.EQ.INX1.AND.DX.LT.0.0) X0=XM
				    CALL PLOT(X0,Y1,3)
				    CALL PLOT(X1,Y1,2)
130				 CONTINUE
			    ELSE
				 CALL PLOT(X,Y1,3)
				 CALL PLOT(XM,Y1,2)
			    ENDIF
150			CONTINUE
		ENDIF
	Y0=IY*DY+Y
	CALL PLOT(X,Y0,3)
	CALL PLOT(XM,Y0,2)
200	CONTINUE
	IF (IF.EQ.1) CALL NEWPEN(0)		! SOLID LINES
	RETURN
	END
C
	SUBROUTINE LINE(X,Y,N,K,J,L,IX,IY,XMIN,DX,YMIN,DY)
C
C	COPYRIGHT (C) CERRITOS COMPUTER SERVICES, INC. 1976
C
C	MODIFIED D.LONG 24-MAY-1983	AT JPL
C
C	X	(R) ARRAY OF UNSCALED ABCISSA VALUES
C	Y	(R) ARRAY OF UNSCALED ORDINATE VALUES
C	N	(I) NUMBER OF POINTS IN THE ARRAY
C	K	(I) REPEAT CYCLE OF A MIXED ARRAY (NORMALLY =1)
C	J	(I) >0, SYMBOL AND LINE SYMBOL EVERY JTH POINT
C		    =0, LINE ONLY
C		    <0, SYMBOL ONLY SYMBOL EVERY JTH POINT
C	L	(I) NUMBER OF SYMBOL ,SEE SYMBOL ROUTINE FOR LIST
C	IX,IY	(I) FIRST REVALENT VALUE IN X,Y ARRAY
C	XMIN,YMIN (R) MINIMUM VALUES IN ARRAY (DETERMINED BY SCALE)
C	DX,DY	(R) SCALED SMOOTHED INCREMENTS (DETERMINED BY SCALE)
C		    IF DX=0.0, THE PLOTTED X IS ZERO
C		    IF DY=0.0, THE PLOTTED Y IS ZERO
C
C		XPLOTTED=(X(I)-XMIN)/DX
C		YPLOTTED=(Y(I)-YMIN)/DY
C
	REAL X(1),Y(1)
	IF (N.LT.1) RETURN
	JA=IABS(J)
	NC=-1
	IF (J.GT.0) NC=-2
	AL=L+.5
	IF (AL.LE.0.0.OR.AL.GT.117.) AL=100.
	I2=1
	I3=3
	IF (J.GE.0) I3=2
	I4=3
	I5=N
100	IF (DX.NE.0.0) THEN
		X1=(X(I2)-XMIN)/DX
	ELSE
		X1=0.0
	ENDIF
	IF (DY.NE.0.0) THEN
		Y1=(Y(I2)-YMIN)/DY
	ELSE
		Y1=0.0
	ENDIF
	CALL PLOT(X1,Y1,I4)
	I4=I3
	IF (J.NE.0) THEN
		IF(MOD(I2-1,JA).EQ.0) CALL 
     $			SYMBOL(X1,Y1,0.14,AL,0.0,NC,-1)
	ELSE
		CALL PLOT(X1,Y1,I4)
	ENDIF
	I2=I2+K
	I5=I5-1
	IF (I5.GT.0) GOTO 100
110	RETURN
	END
C
C
	SUBROUTINE DASHL(X,Y,N,K,J,L,IX,IY,XMIN,DX,YMIN,DY)
C
C	MODIFIED FROM LINE D.LONG 24-MAY-1983 AT JPL
C
C	X	(R) ARRAY OF UNSCALED ABCISSA VALUES
C	Y	(R) ARRAY OF UNSCALED ORDINATE VALUES
C	N	(I) NUMBER OF POINTS IN THE ARRAY
C	K	(I) REPEAT CYCLE OF A MIXED ARRAY (NORMALLY =1)
C	J	(I) >0, SYMBOL AND LINE SYMBOL PLOTTED EVERY JTH POINT
C		    =0, LINE ONLY
C		    <0, SYMBOL ONLY SYMBOL PLOTTED EVERY JTH POINT
C	L	(I) NUMBER OF SYMBOL ,SEE SYMBOL ROUTINE FOR LIST
C	IX,IY	(I) FIRST REVALENT VALUE IN X,Y ARRAY
C	XMIN,YMIN (R) MINIMUM VALUES IN ARRAY (DETERMINED BY SCALE)
C	DX,DY	(R) SCALED (SMOOTHED) INCREMENT (DETERMINED BY SCALE)
C
	REAL X(1),Y(1)
	LOGICAL PEN
	DATA DL/.3/,SL/.15/		! DASH AND SPACE LENGTH
	IF(N.LT.1)GOTO 110
	AL=L+.5
	IF(AL.LE.0.0.OR.AL.GT.117.)AL=100.
	I2=IX
	I2Y=IY
	I3=3
	IF(J.GE.0)I3=2
	IF (DX.NE.0.0) THEN
		X0=(X(I2)-XMIN)/DX
	ELSE
		X0=0.0
	ENDIF
	IF (DY.NE.0.0) THEN
		Y0=(Y(I2Y)-YMIN)/DY
	ELSE
		Y0=0.0
	ENDIF
	CALL PLOT(X0,Y0,3)
	I6=I3
	I4=I3
	I5=N-1
	D=0.0
	PEN=.TRUE.
	IF( J.NE.0) THEN
		IF (MOD((I2-1),J).EQ.0) CALL
     $			SYMBOL(X0,Y0,0.14,AL,0.0,-1,-1)
	ENDIF
	IF (I5.LE.0) GOTO 110
100	IF (DX.NE.0.0) THEN
		X1=(X(I2)-XMIN)/DX
	ELSE
		X1=0.0
	ENDIF
	IF (DY.NE.0.0) THEN
		Y1=(Y(I2Y)-YMIN)/DY
	ELSE
		Y1=0.0
	ENDIF
	D0=((X1-X0)**2+(Y1-Y0)**2)**.5
	I4=3
	DC=SL
	IF (PEN) DC=DL
	IF (D+D0.LT.DC) THEN
		IF (PEN) I4=I3
		CALL PLOT(X1,Y1,I4)
		IF(J.NE.0)THEN
	IF (MOD((I2-1),J).EQ.0)CALL SYMBOL(X1,Y1,0.14,AL,0.0,-1,-1)
		ENDIF
		D=D0+D
		X0=X1
		Y0=Y1
		GOTO 105
	ENDIF
	DD=DC-D
	D=0.0
	A=ATAN2(Y1-Y0,X1-X0)
	X0=X0+DD*COS(A)
	Y0=Y0+DD*SIN(A)
	IF (PEN) I4=I3
	CALL PLOT(X0,Y0,I4)
	PEN=.NOT.PEN
	GOTO 100
105	I2=I2+K
	I2Y=I2Y+K
	I5=I5-1
	IF(I5.GT.0)GOTO 100
110	RETURN
	END
C
C
	SUBROUTINE PLRLN(R,T,N,J,L,IR,RMIN,DR)
C
C	CREATED BY D. LONG 8-AUG-83	AT JPL
C
C	DRAWS A CONNECTED AND/OR SYMBOLED LINE WHEN POINTS ARE SPECIFIED
C	IN POLAR COORDINATES.
C
C	R	(R) ARRAY OF UNSCALED RADIUS VALUES
C	T	(R) ARRAY OF UNSCALED ANGLE VALUES IN DEGREES
C	N	(I) NUMBER OF POINTS IN THE ARRAY
C	J	(I) >0, SYMBOL AND LINE.  SYMBOLS PLOTTED EVERY JTH POINT
C		    =0, LINE ONLY
C		    <0, SYMBOL ONLY. SYMBOLS PLOTTED EVERY JTH POINT
C	L	(I) NUMBER OF SYMBOL ,SEE SYMBOL ROUTINE FOR LIST
C	IR	(I) FIRST REVALENT VALUE IN ARRAYS
C	RMIN	(R) MINIMUM VALUES IN R ARRAY (DETERMINED BY SCALE)
C	DR	(R) SCALED (SMOOTHED) INCREMENT (DETERMINED BY SCALE)
C		    IF DR=0.0, THEN PLOTTED R IS ZERO
C
	REAL R(1),T(1)
	K=1
	TPI=3.141592654/180.
	IF (N.LT.1) GOTO 110
	AL=L+.5
	IF (AL.LE.0.0.OR.AL.GT.117.) AL=100.
	I2=IR
	I3=3
	IF(J.GE.0)I3=2
	I4=3
	I5=N-IR+1
100	IF (DR.NE.0.0) THEN
		R1=(R(I2)-RMIN)/DR
	ELSE
		R1=0.0
	ENDIF
	T1=T(I2)*TPI
	X1=R1*COS(T1)
	Y1=R1*SIN(T1)
	CALL PLOT(X1,Y1,I4)
	I4=I3
	IF(J.NE.0)THEN
	IF (MOD((I2-IX),J*K).EQ.0)CALL SYMBOL(X1,Y1,0.14,AL,0.0,-1,-1)
	ENDIF
	I2=I2+K
	I5=I5-K
	IF(I5.GT.0)GOTO 100
110	RETURN
	END
C
	SUBROUTINE PLRAX(X0,Y0,R0,ASTAR,AEND,AL0,AL1)
C
C	CREATED D. LONG 8-AUG-83	AT JPL
C
C	PLOTS A POLAR COORDINATE SYSTEM AXIS
C
C	X0,Y0	(R) X,Y  COORDINATE VALUES FOR CIRCULAR AXIS CENTER
C	R0	(R) RADIUS OF AXIS
C	ASTAR	(R) STARTING ANGLE OF AXIS, IN DEGREES FROM HORIZONTAL
C	AEND	(R) ENDING ANGLE OF AXIS IN DEGREES FROM HORIZONTAL
C			0.,360. YIELDS FULL CIRCLE PLOT
C	AL0	(R) STARTING ANGLE LABEL NUMERIC VALUE
C	AL1	(R) ENDING ANGLE LABEL NUMERIC VALUE
C	NOTE:	IF AL0=AL1 THEN NO LABEL NUMBERS ARE PLOTTED
C
	TPI=3.141592654/180.
	CS=.15
	XD=CS*2.51
	YD=CS/2.
	N1=R0+.5
	IF (N1.LT.1) RETURN
	RT=FLOAT(N1)+.1
	RD=AEND-ASTAR
	ND=RD/90.+.49
	IF (ASTAR.EQ.AEND) RETURN
	T=ASTAR*TPI
	X1=COS(T)*RT
	Y1=SIN(T)*RT
	CALL PLOT(X0,Y0,3)
	CALL PLOT(X1,Y1,2)
	X2=X1+COS(T)*XD-XD+.05
	Y2=Y1+SIN(T)*(YD+.15)-YD
	IF (AL0.NE.AL1) CALL NUMBER(X2,Y2,CS,AL0,0.,1,-1)
	AC=(AL1-AL0)/(AEND-ASTAR)
	IF (ND*90..GT.AEND) GOTO 150
	DO 100 I=1,ND
		T=(ASTAR+FLOAT(I-1)*90.+45.)*TPI
		CT=COS(T)
		ST=SIN(T)
		DO 50 J=1,N1
			RT1=J-.1
			RT2=J+.1
			X1=CT*RT1
			Y1=ST*RT1
			CALL PLOT(X1,Y1,3)
			X1=CT*RT2
			Y1=ST*RT2
			CALL PLOT(X1,Y1,2)
50		CONTINUE
		IF (ASTAR+360..EQ.AEND.AND.I.EQ.ND) GOTO 175
		T=(ASTAR+FLOAT(I)*90.)*TPI
		CT=COS(T)
		ST=SIN(T)
		X1=CT*RT
		Y1=ST*RT
		CALL PLOT(X1,Y1,3)
		CALL PLOT(X0,Y0,2)
		X2=X1+CT*XD-XD+.05
		Y2=Y1+ST*(YD+.15)-YD
		AD=AC*FLOAT(I)*90.+AL0
		IF (AL1.NE.AL0) CALL NUMBER(X2,Y2,CS,AD,0.,1,-1)
100	CONTINUE
150	IF (ASTAR+FLOAT(ND)*90..NE.AEND) THEN
		T=AEND*TPI
		CT=COS(T)
		ST=SIN(T)
		X1=CT*RT
		Y1=ST*RT
		CALL PLOT(X1,Y1,3)
		CALL PLOT(X0,Y0,2)
		X2=X1+CT*XD-XD+.05
		Y2=Y1+ST*(YD+.15)-YD
		AD=AC*FLOAT(I)*90.+AL0
		IF (AL1.NE.AL0) CALL NUMBER(X2,Y2,CS,AD,0.,1,-1)
	ENDIF
175	DO 200 I=1,N1
		R1=I
		CALL CIRCLE(X0,Y0,ASTAR,AEND,R1,R1,0)
200	CONTINUE
	RETURN
	END
C
	SUBROUTINE CIRCLE(X,Y,A1,A2,R1,R2,D)
C
C	CREATED BY D. LONG    AUG, 1983		AT JPL
C
C	THIS ROUTINE PLOTS A CIRCLE USING THE PLOTTING RESOLUTION
C	AND CIRCLE SIZE TO MAKE A GOOD "APPERANCE" OF A CONTINUOUS CIRCLE.
C
C	X,Y	(R) LOCATION OF CIRCLE CENTER
C	A1	(R) ANGLE OF STARTING POINT FROM HORIZONTAL IN DEGREES
C	A2	(R) ANGLE TO ENDING POINT IN DEGREES
C	R1	(R) RADIUS OF STARTING POINT
C	R2	(R) RADIUS OF ENDING POINT
C	D	(R) D=0.0	SOLID LINE
C		    D=0.5	DASHED LINE
C
	LOGICAL FLAG,PEN
	FLAG=.TRUE.
	TPI=3.141592654/180.0
	SIZE=.1 			! CIRCLE RESOLUTION
	DASHL=.3			! DASH LINE DASH SIZE
	SPACE=DASHL/2.1 		! DASH LINE SPACE SIZE
	A=A1
	R=R1
	ARG=A*TPI
	X1=R*COS(ARG)+X
	Y1=R*SIN(ARG)+Y
	X2=X1
	Y2=Y1
	CALL PLOT(X1,Y1,3)		! PLOT START POINT
	PEN=.TRUE.
	IPEN=2
	DIS=0.
	IF (A1.EQ.A2.OR.R1.LT.0..OR.R2.LT.0.) RETURN
	RD=(R2-R1)/(A2-A1)
100	ARG=A*TPI
	IF (ABS(ACOS(SIZE/R)).LT.ABS(ASIN(SIZE/R))) THEN ! COMPUTE ANGLE INCREM.
		AI=ABS(ACOS(SIZE/R)/TPI)+.01		! BASED ON RESOLUTION
	ELSE						! OF CIRCLE
		AI=ABS(ASIN(SIZE/R)/TPI)+.01
	ENDIF
	A=A+AI*SIGN(1.,A2-A1)
	IF (ABS(A-A2).LT.AI+.01) THEN
		IF (FLAG) THEN				! FINISH LAST POINT
			A=A2
			FLAG=.FALSE.
			IPEN=2
			GOTO 110
		ENDIF
		CALL PLOT(X,Y,3)			! END CURVE
		RETURN
	ENDIF
110	R=R1+RD*(A-A1)
	ARG=A*TPI
	X1=R*COS(ARG)+X
	Y1=R*SIN(ARG)+Y
	CALL PLOT(X1,Y1,IPEN)
	IF (D.EQ.0.5) THEN				! DASHED CURVE
		DIS=SQRT((X1-X2)**2+(Y1-Y2)**2)+DIS
		IF (DIS.GT.DASHL.AND.PEN) THEN
			PEN=.FALSE.
			IPEN=3
			DIS=0.
			GOTO 120
		ENDIF
		IF (DIS.GT.SPACE.AND..NOT.PEN) THEN
			PEN=.TRUE.
			IPEN=2
			DIS=0.
		ENDIF
120		X2=X1
		Y2=Y1
	ENDIF
	GOTO 100
	END
C
	SUBROUTINE LGAXS(X0,Y0,A0,N0,S0,T0,C0,D0,ICOL)
C
C	THIS ROUTINE PLOTS A LOG AXIS WITH EXPONENT LABELING
C	MODIFIED FROM AXIS D.LONG 22-OCT-83	AT JPL
C
C	X0	X COORDINATE OF START OF AXIS
C	Y0	Y COORDINATE OF START OF AXIS
C	A0	CHARACTER STRING TO DESCRIBE AXIS
C	N0	NUMBER OF CHARACTERS IN STRING
C		- ON CLOCKWISE SIDE OF AXIS (NORMAL FOR X)
C		+ ON COUNTER CLOCKWISE SIDE OF AXIS (NORMAL FOR Y)
C		HUNDREDS DIGIT	   = 1 TICKS AND LINE ONLY--NO LABELING
C		THOUSANDS DIGIT    = 1 HORIZONTAL NUMBERS IN Y AXIS LABEL
C		100 THOUSANDS DIGIT = 1 USE COLOR ARRAY
C	S0	LENGTH OF AXIS
C		< 0	TICKS ON OPPOSITE SIDE OF AXIS FROM LABEL
C		= 0	NO ACTION (DUMMY CALL)
C		> 0	TICKS ON SAME SIDE AS AXIS LABEL (NORMAL)
C	T0	ANGLE OF AXIS TO X AXIS OF PAPER
C		0.0 FOR X-AXIS
C		90.0 FOR Y-AXIS
C	C0	COORDINATE OF MINIMUM TICK ON AXIS
C	D0	SCALING FACTOR EXPONENT OF THE FORM (NMAX-NMIN)/LENGTH
C	ICOL	COLOR ARRAY (NEEDED ONLY IF MAG(NO)>=10000)
C		ICOL(1) AXIS COLOR
C		ICOL(2) NUMBERS COLOR
C		ICOL(3) TITLE COLOR (RETURN)
C
	REAL A0(1)
	INTEGER ICOL(3)
	LOGICAL LABELS,COLOR
	DATA CS/.15/		! CHARACTER SIZE
	N1=IABS(N0)
	LABELS=.TRUE.
	HOR=T0
	IF (S0.EQ.0.0) RETURN
	T5=0.1			! TICK LENGTH
	B7=T5+.08		! NUMBER DISTANCE FROM AXIS
	B6=B7
	B8=0.0
	IF(N0.LT.0)GOTO 10
	B3=CS*4.
	B4=CS+0.08
	T2=T0
	GOTO 20
10	B3=-CS*4.
	B4=-T5-CS-.05
	T2=T0
	T5=-T5
20	CONTINUE
	COLOR=.FALSE.
	IF (N1.GE.100000) THEN
		N1=MOD(N1,100000)
		COLOR=.TRUE.
	ENDIF
	IF (N1.GE.10000) N1=MOD(N1,10000)
	IF (N1.GE.1000)THEN
		N1=MOD(N1,1000)
		HOR=0.0
		B4=(ABS(T5)+.05)*SIGN(1.,FLOAT(N0))
		IF (N0.GT.0) B4=B4+3.5*CS
		B6=.5*CS
		B8=(.5*CS+ABS(T5))*SIGN(1.,FLOAT(N0))
		IF (N0.LT.0) B8=B8-CS*1.6
	ENDIF
	IF (N1.GE.100) THEN
		N1=MOD(N1,100)
		LABELS=.FALSE.
	ENDIF
	N2=ABS(S0*D0)+.5				! NUMBER OF TICKS
	XL=ABS(S0)/FLOAT(N2)				! DISTANCE BETWEEN TICKS
	T1=T0*0.017453294
	T3=COS(T1)
	T4=SIN(T1)
	HT3=COS(HOR*.017453294)
	HT4=SIN(HOR*.017453294)
	T6=T5*T3
	T5=T5*T4
C
	IF (COLOR) CALL PLOT(FLOAT(ICOL(1)),0.,0)	!COLOR
	S5=SIGN(1.,S0)*T5
	S6=SIGN(1.,S0)*T6
	X1=X0
	Y1=Y0
	CALL PLOT(X1,Y1,3)
	DO 100 I=1,N2			! PLOT MAJOR TICKS
	DO 100 J=2,10			! PLOT MINOR TICKS
	AT=.6
	IF (J.EQ.2) AT=1.
	X2=X1-S5*AT
	Y2=Y1+S6*AT
	CALL PLOT(X2,Y2,3)
	CALL PLOT(X1,Y1,2)
	AJ=ALOG10(FLOAT(J))-ALOG10(FLOAT(J-1))
	X1=X1+T3*XL*AJ
	Y1=Y1+T4*XL*AJ
	IF (T0.EQ.90.) X1=X0
	CALL PLOT(X1,Y1,2)
100	CONTINUE
	X2=X1-S5
	Y2=Y1+S6
	CALL PLOT(X2,Y2,2)
C
	IF (COLOR) CALL PLOT(FLOAT(ICOL(3)),0.,0)	!COLOR
	IF (.NOT.LABELS) RETURN
	C1=AINT(C0)
	X2=X0-B4*T4-B7*T3		! LOCATE CENTER NUMBER LABELS
	Y2=Y0+B4*T3-B6*T4
	N2=N2+1
	IF (COLOR) CALL PLOT(FLOAT(ICOL(2)),0.,0)	!COLOR
	DO 150 I=1,N2
		CALL SYMBOL(X2,Y2,CS,'10',HOR,2,-1)
		X1=X2-.6*CS*HT4+CS*2.*HT3
		Y1=Y2+.6*CS*HT3+CS*2.*HT4
		IF (C1.EQ.0.0) THEN
		    CALL SYMBOL(X1,Y1,CS*.8,'0',HOR,1,-1)
		ELSE
		    CALL NUMBER(X1,Y1,CS*.8,C1,HOR,-1,-1)
		ENDIF
		C1=C1+SIGN(1.,D0)
		X2=X2+T3*XL
		Y2=Y2+T4*XL
150	CONTINUE
	IF (N1.LT.1) GOTO 155
	X1=0.0
	Y1=0.0
	CALL SYMBOL(X1,Y1,CS,A0,T2,N1,-3)	! GET TITLE LENGTH
	B1=0.5*(ABS(S0)-X1)
	X2=X0+B1*T3-B3*T4-B8*T4
	Y2=Y0+B1*T4+B3*T3
	IF (COLOR) CALL PLOT(FLOAT(ICOL(3)),0.,0)	!COLOR
	CALL SYMBOL(X2,Y2,CS,A0,T2,N1,-1)
155	RETURN
	END
C
C
	SUBROUTINE SCALG(X,S,N,K,IX,XMIN,DX)
C
C	COMPUTES SCALE FACTORS FOR INPUT DATA USING LOG SCALING
C	MODIFIED FROM SCALE D.LONG 4-AUG-83	AT JPL
C
C	X	ARRAY OF DATA TO BE SCANNED FOR MAXIMUM AND MINIMUM
C		VALUES.
C	S	LENGTH OVER WHICH THIS DATA IS TO BE PLOTTED.
C	N	NUMBER OF DATA POINTS IN THE ARRAY X.
C	K	REPEAT CYCLE OF MIXED ARRAY(NORMALLY 1).
C	IX	FIRST RELEVANT DATA POINT IN X
C	XMIN	SMOOTHED MINIMUM AFTER CALL
C	DX	SMOOTHED INCREMENT AFTER CALL
C
C	NOTE: IF X(I)=0.0 THEN 1.E-38 IS USED INSTEAD.
C	      ABSOLUTE VALUE OF X ARRAY IS USED.
C
	REAL X(1)
	NP=N*K
	XI=ABS(X(1))
	IF (XI.EQ.0.0) XI=1.E-38
	XMAX=ALOG10(XI)
	XMIN=XMAX
	DO 100 I=1,NP,K
		XI=ABS(X(I))
		IF (XI.EQ.0.0) XI=1.E-38
		XI=ALOG10(XI)
		XMAX=AMAX1(XMAX,XI)
		XMIN=AMIN1(XMIN,XI)
100	CONTINUE
	SJ=AINT(XMIN)
	IF (XMIN.LT.SJ) SJ=SJ-1.
	SJ2=ANINT(XMAX+.5)
	DX=(SJ2-SJ)/S
	XMIN=SJ
	RETURN
	END
C
	SUBROUTINE LGLIN(X,Y,N,K,J,L,LG,IX,IY,XMIN,DX,YMIN,DY)
C
C	PLOTS A LOGRITHMIC LINE
C	D.LONG 8-AUG-1983	AT JPL
C
C	X	(R) ARRAY OF ABCISSA VALUES
C	Y	(R) ARRAY OF ORDINATE VALUES
C	N	(I) NUMBER OF POINTS IN THE ARRAY
C	K	(I) REPEAT CYCLE OF A MIXED ARRAY (NORMALLY =1)
C	J	(I) >0, SYMBOL AND LINE.  SYMBOLS PLOTTED EVERY JTH POINT
C		    =0, LINE ONLY
C		    <0, SYMBOL ONLY. SYMBOLS PLOTTED EVERY JTH POINT
C	L	(I) NUMBER OF SYMBOL
C	LG	(I) = -2 : X AND Y ARE LOGARITHMIC
C		    = -1 : X IS LOGARITHMIC, Y IS LINEAR
C		    =  1 : X IS LINEAR, Y IS LOGARITHMIC
C	IX,IY	(I) FIRST REVALENT VALUE IN X,Y ARRAY
C	XMIN,YMIN (R) MINIMUM VALUES IN ARRAY (DETERMINED BY SCALG)
C	DX,DY	(R) SCALED (SMOOTHED) INCREMENT (DETERMINED BY SCALG)
C
C	NOTE: IF LOG OF ZERO IS GOING TO OCCUR 1.E-38 IS SUBSTITUTED.
C	      ABSOLUTE VALUE OF X,Y ARRAYS ARE USED WHEN LOG IS COMPUTED.
C
	REAL X(1),Y(1)
	IF (N.LT.1) GOTO 110
	AL=L+.5
	IF (AL.LE.0.0.OR.AL.GT.117.) AL=100.
	I2=1
	I3=3
	IF (J.GE.0) I3=2
	I4=3
	I5=N
100	XI=X(I2)
	YI=Y(I2)
	IF (LG.EQ.1) GOTO 101
	IF (XI.EQ.0.0) XI=1.E-38
	XI=ALOG10(ABS(XI))
101	IF (LG.EQ.-1) GOTO 102
	IF (YI.EQ.0.0) YI=1.E-38
	YI=ALOG10(ABS(YI))
102	IF (DX.NE.0.0) THEN
		X1=(XI-XMIN)/DX
	ELSE
		X1=0.0
	ENDIF
	IF (DY.NE.0.0) THEN
		Y1=(YI-YMIN)/DY
	ELSE
		Y1=0.0
	ENDIF
	CALL PLOT(X1,Y1,I4)
	I4=I3
	IF (J.NE.0) THEN
		IF (MOD((I2-IX),J*K).EQ.0) CALL
     $			SYMBOL(X1,Y1,0.14,AL,0.0,-1,-1)
	ENDIF
	I2=I2+K
	I5=I5-1
	IF (I5.GT.0) GOTO 100
110	RETURN
	END
C
	SUBROUTINE RECT(X0,Y0,X1,Y1)
C
C	ROUTINE TO PLOT A RECTANGLE
C	CREATED BY D. LONG JULY, 1983	AT JPL
C
C	X0,Y0	(R) LOWER LEFT HAND CORNER
C	X1,Y1	(R) UPPER RIGHT HAND CORNER
C
C	NOTE: PEN ENDS UP DOWN AT LOWER-LEFT HAND CORNER
C
	CALL PLOT(X0,Y0,3)
	CALL PLOT(X1,Y0,2)
	CALL PLOT(X1,Y1,2)
	CALL PLOT(X0,Y1,2)
	CALL PLOT(X0,Y0,2)
	RETURN
	END
C
C
	SUBROUTINE SCALE(X,S,N,K,IX,XMIN,DX)
C
C	COPYRIGHT (C) CERRITOS COMPUTER SERVICES, INC. 1976
C
C	CREATES SMOOTHED LINEAR SCALE FACTORS FROM INPUT DATA
C	MODIFIED D.LONG 24-MAY-83	AT JPL
C
C	X	(R) ARRAY OF DATA TO BE SCANNED FOR MAXIMUM AND MINIMUM
C		    VALUES. 
C	S	(R) LENGTH OVER WHICH THIS DATA IS TO BE PLOTTED.
C	N	(I) NUMBER OF DATA POINTS IN THE ARRAY X.
C	K	(I) REPEAT CYCLE OF MIXED ARRAY(NORMALLY 1).
C	IX	(I) FIRST RELEVANT DATA POINT IN X
C	XMIN	(R) SMOOTHED MINIMUM AFTER CALL
C	DX	(R) SMOOTHED INCREMENT AFTER CALL
C
C	TO USE SMOOTHED VALUES: XPLOTTED=(XVALUE-XM)/DX
C
	REAL X(1),Q(6)
	DATA Q/1.0,2.0,4.0,5.0,8.0,10.0/
	NP=N*K
	XMAX=X(1)
	XMIN=XMAX
	DO 100 I=IX,NP,K
		XI=X(I)
		XMAX=AMAX1(XMAX,XI)
		XMIN=AMIN1(XMIN,XI)
100	CONTINUE
	IF (S.LE.0.0) GOTO 160
	DX=(XMAX-XMIN)/S
	IF (DX.LE.0.0) GOTO 160
	SJ=0.0
	IF (DX.LT.1.0) SJ=-1.0
	IDX=ALOG10(DX)+SJ
	DX=DX/(10.0**IDX)
	DO 110 I=1,6
		XI=Q(I)
		IF (XI.GE.DX) GOTO 120
110	CONTINUE
120	DX=XI*(10.0**IDX)
	SI=1.0
	SJ=0.0
	IF (XMIN) 130,170,140
130	SI=-1.0
	SJ=-0.99999
	XMIN=-XMIN
140	IDX=ALOG10(XMIN)+SJ
	XMIN=XMIN/(10.0**IDX)
	XMIN=XMIN-SJ
	XMIN=IFIX(XMIN)*SI*(10.0**IDX)
	GOTO 170
160	DX=1.0
	XMIN=XMIN-0.5
170	RETURN
	END
C
C
	SUBROUTINE SYMBOL(XA,YA,HI,ASC,TH,NC,IPF)
C
C 	EXTENSIVELY REVISED: DGL  AUG, 1987
C	IMPROVED LOWER CASE FONT, ADDED TEXT CONTINUE, ETC.
C	IF ASCII NULL ENCOUNTERED IN INPUT STRING BEYOND FIRST POSITION
C	ROUTINE TERMINATES.
C
C XA 	(R)	STARTING X-COORD. IF XA>998., END OF STRING OF LAST CALL USED
C YA 	(R)	STARTING Y-COORD. IF YA>998., END OF STRING OF LAST CALL USED
C HI 	(R)	HEIGHT OF CHARACTERS IN STRING
C ASC 	(A)	ASCII STRING (BYTE ARRAY)
C TH 	(R)	ANGLE OF ORIENTATION IN DEGREES
C NC	(I)	NUMBER OF CHARACTERS
C		 > 0 NORMAL
C		 -1 PLOT ONLY SINGLE SYMBOL
C		 -2 PLOT LINE THEN SINGLE SYMBOL
C IPF 	(I)	LOCATION CONTROL.  STRING IS PLOTTED WITH (XA,YA) AT:
C		 -3 LOWER LEFT HAND CORNER, END POINTS OF STRING
C		    ARE RETURNED IN XA,YA, NO PLOTTING DONE
C		 -2 LOWER LEFT HAND CORNER, END POINTS OF STRING ARE
C		    RETURNED IN XA,YA.  STRING IS PLOTTED
C		 -1 LOWER LEFT HAND CORNER.  STRING PLOTTED
C		  0 CENTER
C		  1 LOWER RIGHT HAND CORNER
C		  2 NO PLOTTING
C
C	NOTE THAT CODING FOR THE 16 BIT VALUES ASSUMES 2'S COMPLEMENT
C	REPRESENTATION OF INTEGERS.
C
      REAL ASC(1)
      INTEGER*2 SYM(128),STB(416)
      INTEGER*2 SYN1(8) ,SYN9(8), SYN17(8),SYN25(8),SYN33(8),
     $ SYN41(8),SYN49(8),SYN57(8),SYN65(8),SYN73(8),SYN81(8),
     $ SYN89(8),SYN97(8),SYN105(8),SYN113(8),SYN121(8)
      INTEGER*2 STC1(8) ,STC9(8), STC17(8),STC25(8),STC33(8),
     $ STC41(8),STC49(8),STC57(8),STC65(8),STC73(8),STC81(8),
     $ STC89(8),STC97(8),STC105(8),STC113(8),STC121(8),STC129(8),
     $ STC137(8),STC145(8),STC153(8),STC161(8),STC169(8),
     $ STC177(8),STC185(8),STC193(8),STC201(8),STC209(8)
      INTEGER*2 STC217(8),STC225(8),STC233(8),STC241(8),
     $ STC249(8),STC257(8),STC265(8),STC273(8),STC281(8),
     $ STC289(8),STC297(8),STC305(8),STC313(8),STC321(8),
     $ STC329(8),STC337(8),STC345(8),STC353(8),STC361(8),
     $ STC369(8),STC377(8),STC385(8),STC393(8),STC401(8),
     $ STC409(8)
      EQUIVALENCE (SYM(  1),SYN  1(1)),(SYM(  9),SYN  9(1))
      EQUIVALENCE (SYM( 17),SYN 17(1)),(SYM( 25),SYN 25(1))
      EQUIVALENCE (SYM( 33),SYN 33(1)),(SYM( 41),SYN 41(1))
      EQUIVALENCE (SYM( 49),SYN 49(1)),(SYM( 57),SYN 57(1))
      EQUIVALENCE (SYM( 65),SYN 65(1)),(SYM( 73),SYN 73(1))
      EQUIVALENCE (SYM( 81),SYN 81(1)),(SYM( 89),SYN 89(1))
      EQUIVALENCE (SYM( 97),SYN 97(1)),(SYM(105),SYN105(1))
      EQUIVALENCE (SYM(113),SYN113(1)),(SYM(121),SYN121(1))
      EQUIVALENCE (STB(  1),STC  1(1)),(STB(  9),STC  9(1))
      EQUIVALENCE (STB( 17),STC 17(1)),(STB( 25),STC 25(1))
      EQUIVALENCE (STB( 33),STC 33(1)),(STB( 41),STC 41(1))
      EQUIVALENCE (STB( 49),STC 49(1)),(STB( 57),STC 57(1))
      EQUIVALENCE (STB( 65),STC 65(1)),(STB( 73),STC 73(1))
      EQUIVALENCE (STB( 81),STC 81(1)),(STB( 89),STC 89(1))
      EQUIVALENCE (STB( 97),STC 97(1)),(STB(105),STC105(1))
      EQUIVALENCE (STB(113),STC113(1)),(STB(121),STC121(1))
      EQUIVALENCE (STB(129),STC129(1)),(STB(137),STC137(1))
      EQUIVALENCE (STB(145),STC145(1)),(STB(153),STC153(1))
      EQUIVALENCE (STB(161),STC161(1)),(STB(169),STC169(1))
      EQUIVALENCE (STB(177),STC177(1)),(STB(185),STC185(1))
      EQUIVALENCE (STB(193),STC193(1)),(STB(201),STC201(1))
      EQUIVALENCE (STB(209),STC209(1)),(STB(217),STC217(1))
      EQUIVALENCE (STB(225),STC225(1)),(STB(233),STC233(1))
      EQUIVALENCE (STB(241),STC241(1)),(STB(249),STC249(1))
      EQUIVALENCE (STB(257),STC257(1)),(STB(265),STC265(1))
      EQUIVALENCE (STB(273),STC273(1)),(STB(281),STC281(1))
      EQUIVALENCE (STB(289),STC289(1)),(STB(297),STC297(1))
      EQUIVALENCE (STB(305),STC305(1)),(STB(313),STC313(1))
      EQUIVALENCE (STB(321),STC321(1)),(STB(329),STC329(1))
      EQUIVALENCE (STB(337),STC337(1)),(STB(345),STC345(1))
      EQUIVALENCE (STB(353),STC353(1)),(STB(361),STC361(1))
      EQUIVALENCE (STB(369),STC369(1)),(STB(377),STC377(1))
      EQUIVALENCE (STB(385),STC385(1)),(STB(393),STC393(1))
      EQUIVALENCE (STB(401),STC401(1)),(STB(409),STC409(1))
      DATA SYN	1/     0,     8,    20,    26
     1,    31,	  36,	 43,	50/
      DATA SYN	9/    55,    62,    68,    82
     1,    91,	  98,	101,   104/
      DATA SYN 17/   111,   118,   125,   125
     1,   125,	 125,	125,   125/
      DATA SYN 25/   125,   125,   125,   125
     1,   125,	 125,	125,   125/
      DATA SYN 33/   125,   126,   136,   141
     1,   150,	 161,	174,   186/
      DATA SYN 41/   190,   195,   200,   207
     1,   212,	 215,	218,   221/
      DATA SYN 49/   223,   234,   240,   249
     1,   260,	 266,	276,   287/
      DATA SYN 57/   293,   311,   321,   326
     1,   332,	 336,	341,   345/
      DATA SYN 65/   354,   366,   373,   385
     1,   394,	 402,	409,   414/
      DATA SYN 73/   424,   430,   437,   443
     1,   448,	 452,	457,   463/
      DATA SYN 81/   473,   480,   492,   501
     1,   514,	 519,	526,   530/
      DATA SYN 89/   536,   544,   552,   559
     1,   564,	 567,	571,   577/
      DATA SYN 97/   583,   586,   597,   607
     1,   614,	 625,	635,   641/
      DATA SYN105/   654,   660,   668,   675
     1,   681,	 687,	697,   704/
      DATA SYN113/   714,   724,   735,   740
     1,   748,	 754,	762,   766/
      DATA SYN121/   772,   776,   785,   790
     1,   798,	 803,	811,   816/
      DATA STC	1/ 10842,  2088, 11276, -9686
     1, 10842,	8233,  2320,  5131/
      DATA STC	9/ 11044, -9686, 10842,  3080
     1, -9686,	2666,  7256, 18650/
      DATA STC 17/ 26668, -9716, 10842,  2584
     1, 10780, 23258, 10776,  6684/
      DATA STC 25/ -9718, 11336,  3112, 26842
     1,  2092, 22796, -9701, 10330/
      DATA STC 33/  6764, -9718, 11354,  8547
     1, 24872,	2065,  4945, 21260/
      DATA STC 41/ -9693, 11336,  3176, 10826
     1,  7256, 23258, 11304,  3080/
      DATA STC 49/ -9702, 10826, 22746, -9700
     1,  7256, 11368,  3144, 24794/
      DATA STC 57/ 20516, 18708, -9685, 10842
     1,  9057,	4945,-14630,   577/
      DATA STC 65/  2314, 20993, 12849,-14830
     1, 12641,	9075, 16838, 29489/
      DATA STC 73/ 24579, 21540,-14832,  2888
     1,  6932,	8217, 11305,   626/
      DATA STC 81/ 18630, 26924, 12328, 10545
     1,   836,	3083,-14844,   596/
      DATA STC 89/  2049,  8976, 12843,  6441
     1,-14844, 10849,-14798,  8563/
      DATA STC 97/   785, 29126,  4899,-14847
     1, 10826,	4196,  5216, 19142/
      DATA STC105/ 22570,-14820,  4361, 22726
     1,-14820,	 321, 13510, 18630/
      DATA STC113/ 12584, 11315,   780,  2049
     1,-14804,	 833, 12866,-14807/
      DATA STC121/ 12648, 11315,  2084,  1024
     1, 28870, 11316,  6938,  3092/
      DATA STC129/   259,-14840, 13123,  4120
     1,-14828,	 328,  3075,  8988/
      DATA STC137/ 12320,-14796,  7000,  3092
     1,   259,	8200, 13362, 28870/
      DATA STC145/ 11316,   273, 16838,  3075
     1,  6932,	8217, 12584, 11315/
      DATA STC153/  6948,  4185,   264,   710
     1, 11284, 12595,  8232,  7193/
      DATA STC161/ 16838, 20737,-14831,  2625
     1, 25106,-14814,  6211,-14797/
      DATA STC169/  5200,  8292, 16838, 12572
     1, 17094, 20994, 11300, 12595/
      DATA STC177/-14808,   323, 10248, 13105
     1,  5164,	6674,-14804, 12832/
      DATA STC185/  1060,  5200, 12486, 11315
     1,  6948, 23320,  3092,	 3/
      DATA STC193/ 27846, 12595,  2088,   769
     1,-14836,	3075, 13100, 28976/
      DATA STC201/-14847, 13360,  6235,  1088
     1, 12486, 23092,-14824, 13164/
      DATA STC209/ 10289,   264,  5124,-14829
     1, 22576, 29724,-14844,   833/
      DATA STC217/ 12866, 13169, 18630,   769
     1, 13324, 12486,  6260,-14844/
      DATA STC225/ 16432,-14844,  6704,  1076
     1, 12486,	3176,  1140, 18630/
      DATA STC233/   769, 11276, 12595,  2088
     1, 12486, 11315,  6948,-14824/
      DATA STC241/   328,  3075, 13100, 10289
     1, 21000,-14844, 13104,  9260/
      DATA STC249/  6171,  1114, 18630,   769
     1,  5132,	6427, 10272, 13105/
      DATA STC257/-14804, 13424,   626, 28870
     1,   264,	3075,-14796,   624/
      DATA STC265/-14796,   368,   794,-14796
     1, 11272, 28724,  3112,-14844/
      DATA STC273/ 10352,   538, 11380,-14822
     1, 13424,	2092,  1024, 29382/
      DATA STC281/    48,-15358, 12356,   710
     1, 12338, 22724,  7210,  2666/
      DATA STC289/ 19142, 10776,  7256, 28870
     1,-14821,	6216,  8737,   532/
      DATA STC297/  2049,  1124, 12486,  8784
     1,  7203,	 780,  4098, 25798/
      DATA STC305/  6177,   264,-14844,  1140
     1,  8788,	6177,	264,  5122/
      DATA STC313/ 20678,  7188,  8483,  2072
     1,   769, 17094, 13098,  7001/
      DATA STC321/ 16838,  3075, 25652, 12594
     1,  6184,	4625,-14812, 24624/
      DATA STC329/  7203,-14844,   833,  6722
     1, 27161,-14806,	577,  8971/
      DATA STC337/ 13171, 12486,  2156,  1114
     1, 16838, 16899, 12594,  8390/
      DATA STC345/  8536,   538,  9050,  1052
     1,  8390,	8536,  7203,-14844/
      DATA STC353/  6216,  8993,  3100,   259
     1,-14840, 24624, 13106,  7212/
      DATA STC361/  4627,-14816, 13380, 12900
     1, 10289,	4376,  9234,  8390/
      DATA STC369/  8784,-14812,  3075,  4371
     1,  8472,-14812,  9057,  2674/
      DATA STC377/-14845,  2144,   769, 25612
     1,-14844,	 608,-14812,   352/
      DATA STC385/   786,-14812, 24612,-14844
     1,   833, 13324,  6256,  5137/
      DATA STC393/ 24774,    36,-14844,  2370
     1,  6161, 10529,-15054,  4674/
      DATA STC401/ 12898, 16838,  4618,  8731
     1, 12586, 20677,  6937,-14812/
      DATA STC409/ 13360,     4, 28724,-14844
     1,     0,	   0,	  0,	 0/
C
	LOGICAL OUT
	XA1=XA
	YA1=YA
	X0=0.
	Y0=0.
C
C	CONTINUE STRING OUTPUT FROM END OF LAST CALL TO SYMBOL
C
	IF (XA.GT.998.0) XA1=X
	IF (YA.GT.998.0) YA1=Y
	IPF2=IPF
C
	IF (IPF2) 1,2,3
3	X0=-IABS(NC)*HI
	GOTO 1
2	Y0=-HI/2.
	X0=-IABS(NC)*HI/2.
1	H=TH*0.0174532925
	SI=SIN(H)
	CO=COS(H)
	OUT=.TRUE.
	IF (IPF2.EQ.2) RETURN
	IF (IPF2.EQ.-3) THEN
		OUT=.FALSE.
		XX=X
		YY=Y
	ENDIF
	X=XA1+X0*CO-Y0*SI
	Y=YA1+Y0*CO+X0*SI
	H=HI/7.0
	IF (NC.NE.-2) GOTO 4
	IP=2
	IF (OUT) CALL PLOT(XA1,YA1,IP)
4	IP=3
	IF (OUT) CALL PLOT(X,Y,IP)
	IF (NC) 10,20,30
10	NC0=0
C
C	GET FIRST ASCII OR INTEGER VALUE
C
	JCH=ASC(1)+0.1
	GOTO 50
20	CONTINUE
	IF (NC.EQ.-2.AND.OUT) CALL PLOT(XA1,YA1,2)
	IF (IPF2.LT.-1) THEN
		XA=X
		YA=Y
	ENDIF
	RETURN
30	NC0=NC
	I0=0
40	NC0=NC0-1
	IF (NC0.LT.0) GOTO 20
	I0=1+I0
C
C	PXPCGT GETS THE I0TH CHARACTER FROM THE ASC ARRAY AND 
C	PUTS THE RESULT AS AN ASCII VALUE IN JCH
C
	CALL PXPCGT(ASC,I0,JCH)
C
C	TERMINATE ROUTINE IF ASCII NULL ENCOUNTERED
C
	IF (JCH.EQ.0.AND.I0.GT.1) GOTO 20
50	IF (JCH.LT.0.OR.JCH.GE.128) JCH=32
	I3=SYM(1+JCH)
	IF (JCH.GE.32) GOTO 55
	IP=3
	I1=0
	X=X+H*(-2.*CO+3.*SI)
	Y=Y+H*(-3.*CO-2.*SI)
	XW1=X
	YW1=Y
	XW2=X
	YW2=Y
	GOTO 110
55	IF (JCH.NE.112) GOTO 60
	XW2=X+H*2.*SI
	YW2=Y-H*2.*CO
	IF (OUT) CALL PLOT(XW2,YW2,3)
60	IP=2
	I3=1+I3
C
C	DECODE CHARACTER DATA
C
	I4=(I3-1)/2+1
	IV=STB(I4)
	I4=MOD(I3-1,2)*8
	I1=IBITS(IV,I4,8)
C
	I2=MOD(I1,64)
	IF (I1.GE.64) IP=3
	IY=I2/8
	IX=I2-8*IY
	X1=IX
	Y1=IY
	XW1=X+H*(X1*CO-Y1*SI)
	YW1=Y+H*(Y1*CO+X1*SI)
	XW2=XW1
	YW2=YW1
	IF ((JCH.EQ.103.OR.JCH.EQ.112.OR.JCH.EQ.113.OR.
     $		JCH.EQ.121).AND.(I1.LT.192)) THEN
		XW2=XW2+H*2.*SI
		YW2=YW2-H*2.*CO
	ENDIF
110	IF (OUT) THEN
	  CALL PLOT(XW2,YW2,IP)
C MAKE PERIOD SLIGHTLY LARGER THAN ONE POINT
	  IF (JCH.EQ.46.AND.IP.EQ.2) CALL PLOT(XW2+0.0026,YW2,IP)
	ENDIF
	IF (I1.LT.192) GOTO 60
	X=XW1
	Y=YW1
	GOTO 40
	END
C
C
	SUBROUTINE NUMBER(X,Y,HGHT,Z,T,F0,IPF)
C
C	WRITTEN BY D. LONG    AUG, 1983 AT JPL
C	
C	PLOTS A FLOATING POINT NUMBER USING "ENCODE"
C	SEE THE ROUTINE OLDNUMB FOR A FORM OF NUMBER NOT USING ENCODE
C	NOTE THAT THE TEST FOR F0 BEING REAL/INTEGER IS MACHINE DEPENDENT.
C
C	X,Y	COORDINATES OF STRING
C		  (999.,999.) CONTINUE FROM LAST POINT
C	HGHT	HEIGHT OF THE PLOTTED NUMBER
C	Z	FLOATING POINT NUMBER TO BE PLOTTED
C	T	ORIENTATION ANGLE
C	F0	PLOTTING FORMAT (Fn.j)
C		  n = TOTAL NUMBER OF SPACES TO USE (INCLUDING SIGN AND D.P.)
C			[MAX 18 CHARACTERS WIDE]
C		  j = DIGITS TO RIGHT OF DECIMAL POINT (2 DIGITS)
C		  F0 < 0 SPECIFIES EXPONENTIAL NOTATION
C		  F0 = PLOT INTEGER PORTION ONLY WITH NO D.P. (FREE FORMAT)
C		  F0 > 1000 PLOT INTEGER PORTION IN FIXED FORMAT
C		  IF n=0 THEN J IS NUMBER OF DIGITS TO RIGHT OF DECIMAL POINT
C		  WHEN Z OVERFLOWS THIS FORMAT, SPACE IS FILLED WITH ASTERICKS
C	IPF	NUMBER CENTERING FLAG (SEE SYMBOL)
C		=-3 X,Y ARE LOWER LEFT CORNER, END OF STRING RETURNED IN X,Y
C			BUT NUMBER IS NOT PLOTTED
C		=-2 X,Y ARE LOWER LEFT CORNER, END OF STRING RETURNED IN X,Y
C		=-1 X,Y ARE LOWER LEFT CORNER
C		= 0 X,Y ARE STRING CENTER
C		=+1 X,Y ARE LOWER RIGHT CORNER
C		=+2 NO PLOT OUTPUT
C NOTE:
C	F0 IS SET UP TO HANDLE AN INTEGER ARGUMENT (MACHINE DEPENDENT)
C		(NUMBER OF DIGITS TO RIGHT OF DECIMAL POINT)
C
	BYTE B(18)			! WORKING BUFFER
	EQUIVALENCE (NN,RR)		! FOR REAL/INTEGER TEST
C
	ND=-1
	HG=HGHT
	T1=T
	FA=F0
	IF (FA.NE.0.0) GOTO 3
	RR=FA
	IF (NN.NE.0) GOTO 2		! CHECK FOR INTEGER F0 INPUT
    	ND=0
	FA=1.0
	GOTO 10
  1	NN=AMOD(FA,1000.)		! PLOT FORMATED INTEGER 
	FA=0.
	GOTO 10
  2	FA=1.0
	IF (NN.LT.0) FA=0.0
	ND=NN
	NN=0
	GOTO 10
  3	IF (ABS(FA).GT.1.E30)  THEN
		RR=FA
		GOTO 2			! INTEGER -1 INPUT
	ENDIF
	IF (FA.GT.999) GOTO 1
	F=ABS(FA)*1.000002		! GET PLOT FORMAT
	NN=F
	F=(F-NN)*100.
	ND=F
 10	IF (ND.GT.17) ND=ND/10		! SOLVE SIMPLE ERRORS
	IF (NN.EQ.0) THEN		! DIGITS TO LEFT OF DECIMAL POINT
		NN=ND+2
		IF (Z.EQ.0.AND.FA.EQ.0.0) NN=1
		IF (Z.NE.0.) THEN
			ALG=ALOG10(ABS(Z))
			IF (ALG.LT.0.0) ALG=0.
			NN=ND+2+ALG
		ENDIF
		IF (Z.LT.0.0) NN=NN+1
		IF (FA.LT.0.0) NN=NN+4
	ENDIF
	IF (ND.GT.NN) GOTO  60		! FORMAT ERROR
	IF (NN.GT.18) NN=18		! MAX CHARACTERS
	IF (FA) 20,30,40
 20	ENCODE(18,21,B,ERR=90) Z
 21	FORMAT(E<NN>.<ND>)
	GOTO 50
 30	I=Z
	ENCODE(18,31,B,ERR=90) I
 31	FORMAT(I<NN>)
	GOTO 50
 40	ENCODE(18,41,B,ERR=90) Z
 41	FORMAT(F<NN>.<ND>)
 50	CALL SYMBOL(X,Y,HG,B,T1,NN,IPF)
 60	RETURN
 90	DO 95 I=1,18
		B(I)='*'
		IF (I.EQ.NN-ND) B(I)='.'
 95	CONTINUE
	GOTO 50
	END
C
C
	SUBROUTINE OLDNUMB(X,Y,HGHT,Z,T,N,IPF)
C
C	COPYRIGHT (C) CERRITOS COMPUTER SERVICES, INC. 1976
C	MODIFIED BY DGL    AUG, 1983	AT JPL
C	
C	PLOTS A FLOATING POINT NUMBER.
C	DOES NOT USE "ENCODE" NOR MACHINE DEPENDENT INT/REAL
C	VARIABLE DETECTION.
C
C	X,Y	COORDINATES OF STRING
C	HGHT	HEIGHT OF THE PLOTTED NUMBER IN FLT. PT.
C	Z	FLT. PT. NUMBER TO BE PLOTTED
C	T	ORIENTATION ANGLE FOR THE NUMBER
C	N	NUMBER OF DECIMAL DIGITS FOR OUTPUT.
C		N=-1 WILL SUPPRESS THE DECIMAL POINT.
C	IPF	<0 X,Y ARE LOWER LEFT CORNER
C		=0 X,Y ARE CENTER
C		>0 X,Y ARE LOWER RIGHT CORNER
C
	DATA KZERO/48/
	TPI=3.141592654/180.
	FPN=ABS(Z)*1.000002
	XT=X
	YT=Y
	ACNT=0.
	IPF2=IPF
	IF (IPF2) 1,2,2
2	IF (Z) 4,5,6
5	ACNT=2.+FLOAT(N)
	GOTO 3
4	ACNT=1.
6	IF (ABS(Z)-1.)5,11,11
11	ACNT=AINT(ALOG10(ABS(Z)))+FLOAT(N+1)+ACNT
3	CO=COS(TPI*T)
	SI=SIN(TPI*T)
	X0=0.
	Y0=0.
	IF (IPF2)1,7,8
8	X0=-ACNT*HGHT
	GOTO 9
7	X0=-ACNT*HGHT/2.
	Y0=-HGHT/2.
9	XT=X+X0*CO-Y0*CO
	YT=Y+Y0*SI+X0*SI
1	IF(Z)10,90,20
  10	CALL SYMBOL(XT,YT,HGHT,'-',T,1,-1)
	XT=999.
	YT=999.
  20	I=ALOG10(FPN)+1.0
	IF (I) 50,50,30
  30	DO 40 J = 1,I
	K = FPN*10.0**(J-I)
	CALL SYMBOL(XT,YT,HGHT,FLOAT(K+KZERO),T,-1,-1)
	XT=999.
	YT=999.
	FPN = FPN - FLOAT(K*10**(I-J))
  40	CONTINUE
	IF (N+1)80,80,50
  50	CALL SYMBOL(XT,YT,HGHT,'.',T,1,-1)
	XT=999.
	YT=999.
	IF(N)80,80,60
  60	DO 70 IJ=1,N
	K=FPN*10.0
	CALL SYMBOL(XT,YT,HGHT,FLOAT(K+KZERO),T,-1,-1)
  70	FPN=FPN*10.0-FLOAT(K)
  80	RETURN
  90	CALL SYMBOL(XT,YT,HGHT,'0.0000000000 ',T,2+N,-1)
	RETURN
	END
C
C
	SUBROUTINE PXPCGT(L,I,IV)
C
C	THIS ROUTINE GETS THE ITH BYTE FROM THE ARRAY OF DATA VALUES
C	IN L AND PUTS THE RESULT IN IV.  THIS ROUTINE MAY NEED TO
C	BE MODIFIED FOR DIFFERENT MACHINES.
C
C L	ARRAY HOLDING DESIRED BYTES
C I	POINTER INTO STRING
C IV	VARIABLE TO RECEIVE I0'TH BYTE OF L0
C
	BYTE L(1),L2(4)
	EQUIVALENCE (L2(1),I2)
	I2=0
	L2(1)=L(I)
	IV=I2
	RETURN
	END
C
	SUBROUTINE SHADE(A,B,N,I,L,D,T,W,MA,AM,DA,BM,DB)
C
C	SHADES AN AREA DEFINED BY THE A AND B ARRAYS
C
C	COPYRIGHT (C) CERRITOS COMPUTER SERVICES, INC. 1981
C	MODIFIED AUG 1983 AND 1987 BY D. LONG
C
C A    (R) ARRAY CONTAINING ABSCISSA VALUES
C B    (R) ARRAY CONTAINING ORDINATE VALUES
C N    (I) NUMBER OF POINTS IN ARRAY
C I    (I)INCREMENT BETWEEN POINTS IN ARRAY
C L    (I) SHADE FORMAT CONTROL
C		-3 CLEAR AREA AND OUTLINE
C		-2 CLEAR AREA
C		-1 CLEAR OUTLINE
C		 0 NO ACTION
C		+1 SHADE OUTLINE
C		+2 SHADE AREA
C		+3 SHADE AREA AND OUTLINE
C D    (D) DISTANCE BETWEEN SHADING LINES
C T    (D) ANGLE OF SHADING LINES
C W    (D) WORKING ARRAY DIMENSIONED AT LEAST (3*N)
C MA   (I) TYPE OF SHADING (LINE TYPE) IF MA => 0
C AM   (R) MINIMUM VALUE FOR A ARRAY
C DA   (R) SCALE FACTOR FOR A ARRAY
C BM   (R) MINIMUM VALUE FOR B ARRAY
C DB   (R) SCALE FACTOR FOR B ARRAY
C
	INTEGER KP(2)
	REAL A(1),B(1),W(3,1)
C
C	EQUIVALENCES ARE FOR CONVENIENCE ONLY
C
	EQUIVALENCE (KP(1),KD),(KP(2),KU)
	DATA KP/2,3/
C
	ROT(X,Y)=X*CT+Y*ST
	SCL(Q,R,S)=(Q-R)/S
C
	IF (N.LE.2) GOTO 110
	IF (I.LE.0) GOTO 110
	AO=AM
	AI=DA
	BO=BM
	BI=DB
	IF (AI.LE.0.0) AI=1.0
	IF (BI.LE.0.0) BI=1.0
	LT=IABS(L)
	KD=2
	IF (L.LT.0) KD=9
	IF (LT.LE.1) GOTO 90
	DT=D
	IF (DT.LE.0.0) DT=0.005
	TH=T-(AINT(T/360.0))*360.0
	IF (TH.LT.0.0) TH=TH+360.0
	TH=TH*0.0174532
	ST=SIN(TH)
	CT=COS(TH)
	YS=ROT(SCL(B(1),BO,BI),-SCL(A(1),AO,AI))
	YL=YS
	I1=1
	DO 20 I0=1,N
		TX0=SCL(A(I1),AO,AI)
		TY0=SCL(B(I1),BO,BI)
		W(1,I0)=ROT(TX0,TY0)
		TY1=ROT(TY0,-TX0)
		W(2,I0)=TY1
		YS=AMIN1(YS,TY1)
		YL=AMAX1(YL,TY1)
		I1=I1+I
20	CONTINUE
	YC=YS+0.5*DT
	IF (MA.GE.O) CALL NEWPEN(MA)
30	IF (YC.GT.YL) GOTO 90
	I1=0
	DO 60 I0=1,N
		I2=I0+1
		IF(I2.GT.N)I2=1
		YT0=W(2,I0)
		YT1=W(2,I2)
		W1=AMIN1(YT0,YT1)
		IF(YC.LT.W1)GOTO 60
		W2=AMAX1(YT0,YT1)
		IF (YC.GE.W2) GOTO 60
		XT0=W(1,I0)
		XT1=W(1,I2)
		X1=XT0-(YT0-YC)*(XT1-XT0)/(YT1-YT0)
		I1=I1+1
		I2=I1
40		IF (I2.EQ.1) GOTO 50
		I3=I2-1
		IF (W(3,I3).GE.X1) GOTO 50
		W(3,I2)=W(3,I3)
		I2=I3
		GOTO 40
50		W(3,I2)=X1
60	CONTINUE
	IF (I1.LE.0) GOTO 80
	I2=1
	DO 70 I0=1,I1
		XT0=W(3,I0)
		XT=ROT(XT0,-YC)
		YT=ROT(YC,XT0)
		I2=3-I2
		KC=KP(I2)
		CALL PLOT(XT,YT,KC)
70	CONTINUE
80	YC=YC+DT
	GOTO 30
90	IF (2*(LT/2).EQ.LT) GOTO 110
	CALL PLOT(SCL(A(1),AO,AI),SCL(B(1),BO,BI),KU)
	I1=1
	DO 100 I0=2,N
		I1=I1+I
		CALL PLOT(SCL(A(I1),AO,AI),SCL(B(I1),BO,BI),KD)
100	CONTINUE
	CALL PLOT(SCL(A(1),AO,AI),SCL(B(1),BO,BI),KD)
C
C	RESET PEN TYPE
C
110	IF (MA.GE.O) CALL NEWPEN(0)
	RETURN
	END
C
C
	SUBROUTINE CSHADE(AM,BM,RAD,A1,A2,RES,L,D,T,W,MA1,MA2)
C
C	SHADES A CIRCULAR SEGMENT OR FULL CIRLE
C	WRITTEN AUG 1987 BY D. LONG
C
C AM,BM(R) CENTER OF CIRCLE
C RAD  (R) RADIUS OF CIRCLE SEGMENT
C RES  (R) ALONG ARC RESOLUTION
C A1,A2(R) START,END ANGLES FOR SEGMENT.  IF A2-A1 >= 360, FULL CIRCLE USED
C	    A2 => A1
C L    (I) SHADE FORMAT CONTROL
C		-3 CLEAR AREA AND OUTLINE
C		-2 CLEAR AREA
C		-1 CLEAR OUTLINE
C		 0 NO ACTION
C		+1 DRAW OUTLINE
C		+2 SHADE AREA
C		+3 SHADE AREA AND OUTLINE
C D    (D) DISTANCE BETWEEN SHADING LINES
C T    (D) ANGLE OF SHADING LINES
C W    (D) WORKING ARRAY DIMENSIONED AT LEAST 
C		3*INT((A2-A1)/(180*ATAN(RES/RAD)/PI))+9
C MA1   (I) LINE TYPE OF SHADING (IF MA1 => 0)
C MA2   (I) LINE TYPE OF OUTLINE (IF MA2 => 0)
C
	INTEGER KP(2)
	REAL W(3,1)
	DATA KP/2,3/,DTR/0.0174532/
C
	ROT(X,Y)=X*CT+Y*ST
C
	IF (RAD.LE.0.0) GOTO 110
	AO=AM
	BO=BM
	LT=IABS(L)
	KP(1)=2
	IF (L.LT.0) KP(1)=9
	ANGINC=60.0
	IF (RES.GT.0.0) ANGINC=ATAN(RES/RAD)/DTR
	IF (LT.LE.1) GOTO 90
	DT=D
	IF (DT.LE.0.0) DT=0.005
	TH=T-(AINT(T/360.0))*360.0
	IF (TH.LT.0.0) TH=TH+360.0
	TH=TH*DTR
	ST=SIN(TH)
	CT=COS(TH)
C
C	FILL WORKING ARRAY AND COMPUTE MAX/MIN
C
	I1=1
	ANG=A1
10	CONTINUE
		TX0=RAD*COS(ANG*DTR)+AM
		TY0=RAD*SIN(ANG*DTR)+BM
		W(1,I1)=ROT(TX0,TY0)
		TY1=ROT(TY0,-TX0)
		W(2,I1)=TY1
		IF (I1.EQ.1) THEN
			YS=TY1
			YL=TY1
		ENDIF
		YS=AMIN1(YS,TY1)
		YL=AMAX1(YL,TY1)
		I1=I1+1
		IF (ANG.EQ.A2) GOTO 20
		ANG=ANG+ANGINC
		IF (ANG.GT.A2) ANG=A2
		GOTO 10
20	CONTINUE
C
C	ADD SEGMENT LINES IF REQUIRED
C
	IF (A2-A1.LT.360.0) THEN
		TX0=AM
		TY0=BM
		W(1,I1)=ROT(TX0,TY0)
		TY1=ROT(TY0,-TX0)
		W(2,I1)=TY1
		YS=AMIN1(YS,TY1)
		YL=AMAX1(YL,TY1)
		I1=I1+1
		TX0=RAD*COS(A1*DTR)+AM
		TY0=RAD*SIN(A1*DTR)+BM
		W(1,I1)=ROT(TX0,TY0)
		TY1=ROT(TY0,-TX0)
		W(2,I1)=TY1
		YS=AMIN1(YS,TY1)
		YL=AMAX1(YL,TY1)
		I1=I1+1
	ENDIF
C
C	SHADE AREA
C
	N=I1-1
	YC=YS+0.5*DT
	IF (MA1.GE.0) CALL NEWPEN(MA1)
30	IF (YC.GT.YL) GOTO 90
	I1=0
	DO 60 I0=1,N
		I2=I0+1
		IF (I2.GT.N) I2=1
		YT0=W(2,I0)
		YT1=W(2,I2)
		W1=AMIN1(YT0,YT1)
		IF (YC.LT.W1) GOTO 60
		W2=AMAX1(YT0,YT1)
		IF (YC.GE.W2) GOTO 60
		XT0=W(1,I0)
		XT1=W(1,I2)
		X1=XT0-(YT0-YC)*(XT1-XT0)/(YT1-YT0)
		I1=I1+1
		I2=I1
40		IF (I2.EQ.1) GOTO 50
		I3=I2-1
		IF (W(3,I3).GE.X1) GOTO 50
		W(3,I2)=W(3,I3)
		I2=I3
		GOTO 40
50		W(3,I2)=X1
60	CONTINUE
	IF (I1.LE.0) GOTO 80
	I2=1
	DO 70 I0=1,I1
		XT0=W(3,I0)
		XT=ROT(XT0,-YC)
		YT=ROT(YC,XT0)
		I2=3-I2
		KC=KP(I2)
		CALL PLOT(XT,YT,KC)
70	CONTINUE
80	YC=YC+DT
	GOTO 30
90	IF (2*(LT/2).EQ.LT) GOTO 110
C
C	PLOT AREA OUTLINE
C
	IF (MA2.GE.0) CALL NEWPEN(MA2)
	ANG=A1
	IP=KP(2)
100	CONTINUE
		TX0=RAD*COS(ANG*DTR)+AM
		TY0=RAD*SIN(ANG*DTR)+BM
		CALL PLOT(TX0,TY0,IP)
		IP=KP(1)
		IF (ANG.EQ.A2) GOTO 105
		ANG=ANG+ANGINC
		IF (ANG.GT.A2) ANG=A2
		GOTO 100
105	CONTINUE
C
C	ADD SEGMENT LINES IF REQUIRED
C
	IF (A2-A1.LT.360.0) THEN
		CALL PLOT(AM,BM,IP)
		TX0=RAD*COS(A1*DTR)+AM
		TY0=RAD*SIN(A1*DTR)+BM
		CALL PLOT(TX0,TY0,IP)
	ENDIF
C
C	ALL DONE
C
110	RETURN
	END
C
C
	SUBROUTINE ELLIPSE(X0,Y0,SMAJ,SMIN,ANG,AST,AEND)
C
C	WRITTEN BY:	DG LONG,  JPL  JUNE 1985
C	MODIFIED: DGL NOV 7,1985 TO HANDLE DEGENERATE ELLIPSES
C
C	THIS ROUTINE PLOTS AN ELLIPSE.  THE NUMBER OF LINE SEGMENTS
C	USED TO PLOT THE CURVE IS DEPENDENT ON THE RADIUS OF THE ELLIPSE
C	AT EACH POINT.
C
C	X0,Y0	(R)	CENTER OF ELLIPSE (HALFWAY BETWEEN FOCI)
C	SMAJ	(R)	LENGTH OF SEMI-MAJOR AXIS (DISTANCE FROM ELLIPSE
C				CENTER TO CURVE THROUGH ONE FOCUS)
C	SMIN	(R)	LENGTH OF SEMI-MINOR AXIS (PERPENDICULAR TO SMAJ)
C				0 < SMIN < SMAJ
C	ANG	(R)	ANGLE FROM HORIZONTAL OF SEMI-MAJOR AXIS (DEG)
C	AST	(R)	START ANGLE OF ELLIPSE CURVE FROM SEMI-MAJOR AXIS (DEG)
C	AEND	(R)	END ANGLE OF ELLIPSE CURVE FROM SEMI-MAJOR AXIS (DEG)
C				AEND > AST
C
C	AST AND AEND ARE DEFINED USING THE IDENTITY FOR AN ELLIPSE:
C		(X,Y) = (SMAJ * COS(T) , SMIN * SIN(T) )
C
	DATA ERR/0.025/		! WORST-CASE PLOTTING ERROR
C
	IP=3
	RANG=ANG*3.141592654/180.
	CO=COS(RANG)
	SI=SIN(RANG)
	DTHETA=ASIN(ERR/SMAJ)
	TANG=AST*3.141592654/180.
	EANG=AEND*3.141592654/180.
C
100	CONTINUE
		X=COS(TANG)*SMAJ
		Y=SIN(TANG)*SMIN
		X1=CO*X-SI*Y+X0
		Y1=SI*X+CO*Y+Y0
		CALL PLOT(X1,Y1,IP)
		IP=2
		R=SQRT(X**2+Y**2)
		AA=2.
		IF (R.GT.0.0) AA=ERR/R
		IF (ABS(AA).LT.1.0) THEN
			AA=ASIN(AA)
		ELSE
			AA=0
		ENDIF
		DTHETA=MAX(AA,0.00001)	! MINIMUM STEP SIZE
		TANG=TANG+DTHETA
		IF (TANG.GT.EANG) GOTO 200
	GOTO 100
C
200	CONTINUE
	X=COS(EANG)*SMAJ
	Y=SIN(EANG)*SMIN
	X1=CO*X-SI*Y+X0
	Y1=SI*X+CO*Y+Y0
	CALL PLOT(X1,Y1,2)
	CALL PLOT(X0,Y0,3)	! LEAVE PEN AT ELLIPSE CENTER
C
	RETURN
	END
C
C
	SUBROUTINE ARROW(X,Y,ALEN,ANG,APT,APTAN)
C
C	CREATED BY D. LONG, DEC. 1983 AT JPL
C	PLOTS AN ARROW WITH CHOICE OF HEADS
C
C	X,Y	(R)	LOCATION OF TAIL (INCHES)
C	ALEN	(R)	LENGTH OF ARROW (INCHES)
C	ANG	(R)	ANGLE FROM HORIZONTAL (DEGREES)
C	APT	(R)	LENGTH OF ARROW POINT (INCHES)
C			> 0   -->  (NORMAL ARROW)
C			< 0   -I> (DELTA HEAD ARROW)
C	APTAN	(R)	ANGLE OF ARROW POINT FROM LINE (DEGREES)
C
	DATA API/0.0174532/		! PI/180
	IF (ALEN.EQ.0.0) RETURN
	X0=X+ALEN*COS(API*ANG)
	Y0=Y+ALEN*SIN(API*ANG)
	CALL PLOT(X,Y,3)
	CALL PLOT(X0,Y0,2)
	X1=X0+ABS(APT)*COS(API*(180.-APTAN+ANG))
	Y1=Y0+ABS(APT)*SIN(API*(180.-APTAN+ANG))
	CALL PLOT(X1,Y1,2)
	IF (APT.GT.0.0) CALL PLOT(X0,Y0,2)
	X1=X0+ABS(APT)*COS(API*(APTAN+ANG-180.))
	Y1=Y0+ABS(APT)*SIN(API*(APTAN+ANG-180.))
	CALL PLOT(X1,Y1,2)
	IF (APT.LT.0.0) CALL PLOT(X0,Y0,2)
	CALL PLOT(X,Y,3)
	RETURN
	END
C
C
	SUBROUTINE JPLTAG(X,Y,H,A,IOPT,D,M,SA,W)
C
C	CREATED BY D. LONG  MAR. 1984	AT JPL
C
C	PLOTS THE JPL LOGO
C
C	X,Y	(R)	LOCATION OF LOWER-LEFT CORNER OF JPL
C	H	(R)	HEIGHT (IN INCHES) OF LETTERS
C	A	(R)	ANGLE OF LETTERS
C	IOPT	(I)	OPTION FLAG
C			(ONE'S DIGIT)	= 0 NO ACTION
C					= 1 DRAW OUTLINE
C					= 2 SHADE AREA
C					= 3 SHADE AREA AND OUTLINE
C			(TEN'S DIGITS)	COLOR TABLE VALUE 
C	D	(R)	DISTANCE BETWEEN SHADING LINES
C	M	(R)	LINE TYPE OF SHADING (SEE SHADE)
C	SA	(R)	SHADING ANGLE (DEGREES)
C	W	(R)	WORKING ARRAY DIMENSIONED W(56+54*H)
C
	DIMENSION W(1)			! WORKING ARRAY
C
	ICURVE=3*H+1
	IF (ICURVE.GT.10) ICURVE=10
	AN=3.141592654/FLOAT(2*ICURVE)
	S=ABS(H)/2.2
	IOFF=2*ICURVE+10
	IOPT2=MOD(IOPT,10)
	COLOR=FLOAT(IOPT/10)
	IF (COLOR.NE.0.0) CALL PLOT(COLOR,0.,0)
	SINA=SIN(A*3.141592654/180.)
	COSA=COS(A*3.141592654/180.)
C
C	MAKE J FIRST
C
	W(1)=0.
	W(IOFF+1)=0.
	W(2)=.2
	W(IOFF+2)=.3
	W(3)=.8
	W(IOFF+3)=.3
	DO 1 I=1,ICURVE
		W(3+I)=.2*SIN(AN*I)+.8
		W(3+I+IOFF)=.5-.2*COS(AN*I)
1	CONTINUE
	W(4+ICURVE)=1.0
	W(4+ICURVE+IOFF)=2.2
	W(5+ICURVE)=1.3
	W(5+ICURVE+IOFF)=2.2
	W(6+ICURVE)=1.3
	W(6+ICURVE+IOFF)=.3
	DO 2 I=1,ICURVE
		W(6+I+ICURVE)=1.+.3*COS(AN*I)
		W(6+I+ICURVE+IOFF)=.3-.3*SIN(AN*I)
2	CONTINUE
	INUM=6+2*ICURVE
	DO 10 I=1,INUM
		X1=(W(I)*COSA-W(I+IOFF)*SINA)*S+X
		W(I+IOFF)=(W(I)*SINA+W(I+IOFF)*COSA)*S+Y
		W(I)=X1
10	CONTINUE
	CALL SHADE(W(1),W(1+IOFF),INUM,1,IOPT2,D,SA,W(IOFF+INUM+2),
     $		M,0.,1.,0.,1.)
C
C	MAKE P NEXT
C
	IOFF=4*ICURVE+10
	W(1)=1.7
	W(IOFF+1)=0.
	W(2)=1.7
	W(IOFF+2)=2.2
	W(3)=2.9
	W(IOFF+3)=2.2
	DO 3 I=1,ICURVE
		W(3+I)=2.9+.3*SIN(AN*I)
		W(3+I+IOFF)=1.9+.3*COS(AN*I)
3	CONTINUE
	W(4+ICURVE)=3.2
	W(4+ICURVE+IOFF)=1.3
	DO 4 I=1,ICURVE
		W(4+I+ICURVE)=2.9+.3*COS(AN*I)
		W(4+I+ICURVE+IOFF)=1.3-.3*SIN(AN*I)
4	CONTINUE
	W(5+2*ICURVE)=2.2
	W(5+2*ICURVE+IOFF)=1.0
	W(6+2*ICURVE)=2.0
	W(6+2*ICURVE+IOFF)=1.3
	W(7+2*ICURVE)=2.7
	W(7+2*ICURVE+IOFF)=1.3
	DO 5 I=1,ICURVE
		W(7+I+2*ICURVE)=.2*SIN(AN*I)+2.7
		W(7+I+2*ICURVE+IOFF)=1.5-.2*COS(AN*I)
5	CONTINUE
	W(8+3*ICURVE)=2.9
	W(8+3*ICURVE+IOFF)=1.7
	DO 6 I=1,ICURVE
		W(8+I+3*ICURVE)=2.7+.2*COS(AN*I)
		W(8+I+3*ICURVE+IOFF)=1.7+.2*SIN(AN*I)
6	CONTINUE
	W(9+4*ICURVE)=2.0
	W(9+4*ICURVE+IOFF)=1.9
	W(10+4*ICURVE)=2.0
	W(10+4*ICURVE+IOFF)=0.0
	INUM=10+4*ICURVE
	DO 11 I=1,INUM
		X1=(W(I)*COSA-W(I+IOFF)*SINA)*S+X
		W(I+IOFF)=(W(I)*SINA+W(I+IOFF)*COSA)*S+Y
		W(I)=X1
11	CONTINUE
	CALL SHADE(W(1),W(1+IOFF),INUM,1,IOPT2,D,SA,W(IOFF+INUM+2),
     $		M,0.,1.,0.,1.)
C
C	MAKE L LAST
C
	IOFF=2*ICURVE+8
	W(1)=3.5
	W(IOFF+1)=0.3
	W(2)=3.5
	W(IOFF+2)=2.2
	W(3)=3.8
	W(IOFF+3)=2.2
	W(4)=3.8
	W(IOFF+4)=0.5
	DO 7 I=1,ICURVE
		W(4+I)=4.0-.2*COS(AN*I)
		W(4+I+IOFF)=0.5-.2*SIN(AN*I)
7	CONTINUE
	W(5+ICURVE)=4.6
	W(5+ICURVE+IOFF)=0.3
	W(6+ICURVE)=4.8
	W(6+ICURVE+IOFF)=0.
	W(7+ICURVE)=3.8
	W(7+ICURVE+IOFF)=0.
	DO 8 I=1,ICURVE
		W(7+I+ICURVE)=3.8-.3*SIN(AN*I)
		W(7+I+ICURVE+IOFF)=.3-.3*COS(AN*I)
8	CONTINUE
	INUM=7+2*ICURVE
	DO 14 I=1,INUM
		X1=(W(I)*COSA-W(I+IOFF)*SINA)*S+X
		W(I+IOFF)=(W(I)*SINA+W(I+IOFF)*COSA)*S+Y
		W(I)=X1
14	CONTINUE
	CALL SHADE(W(1),W(1+IOFF),INUM,1,IOPT2,D,SA,W(IOFF+INUM+2),
     $		M,0.,1.,0.,1.)
	RETURN
	END
C
C
	SUBROUTINE PLT3D(A,MD,ND,M,N,D,L,W2,NN,ALT,AZ,XLEN,XOFF,
     $		YLEN,YOFF,ZFAC,ZOFF,IER)
C
C	ROUTINE TO PLOT A GRID SCRIBED ON A 3D SURFACE USING NXTVU
C	TO PROVIDE HIDDEN LINE REMOVAL
C	MODIFIED 8/11/87 BY D.LONG, JPL
C
C	A	(R): MD BY ND ARRAY CONTAINING Z VALUES TO BE PLOTTED
C	MD,ND	(I): DIMENSIONS OF A
C	M,N	(I): SIZE OF DATA IN A ARRAY  M,N > 1
C	D	(R): WORKING ARRAY OF DIMENSIONED D(L)
C	L	(I): DIMENSION OF D ARRAY
C	W2	(R): WORKING ARRAY OF DIMENSIONED W2(NN)
C	NN	(I): DIMENSION OF W2 ARRAY
C	ALT,AZ	(R): VIEWING ALTITUDE, AZIMUTH IN DEGREES
C	XLEN	(R): LENGTH OF UNPROJECTED X AXIS IN INCHES
C			> 0 SURFACE ABOVE PLANE SHOWN
C			< 0 SURFACE BELOW PLANE SHOWN
C	XOFF	(R): ORIGIN OF PLOT
C	YLEN	(R): LENGTH OF UNPROJECTED Y AXIS IN INCHES
C	YOFF	(R): ORIGIN OF PLOT
C	ZFAC	(R): Z COORDINATE SCALE FACTOR (Zplotted=ZFAC*A(I,J)+ZOFF)
C	ZOFF	(R): Z AXIS OFFSET VALUE
C	IER	(I): RETURNED ERROR CODE
C			0: NO ERROR
C			1: WORKING STORAGE W2 INTERNAL SPACE EXCEEDED
C			2: L LESS THAN 2*MIN(M,N)
C
	DIMENSION D(L),A(MD,ND),W2(NN)
C
	COMMON /PLT3B/A1,A2,A3,B1,B2,B3,B4
	DATA DTR/0.0174532925/
C
	IX(I)=(I-1)*2+1
	IY(I)=I*2
C
	LMAX=2*MIN(M,N)
	IF (L.LT.2*LMAX) GOTO 500
	TAZ=AZ*DTR
	TALT=ALT*DTR
	SAZ=SIN(TAZ)
	CAZ=COS(TAZ)
	SAL=SIN(TALT)
	CAL=COS(TALT)
	XSC=ABS(XLEN)/FLOAT(N-1)
	ID=1
	IF (XLEN.LT.0.0) ID=-1
	YSC=YLEN/FLOAT(M-1)
	A1=CAZ*XSC
	A2=-SAZ*YSC
	A3=XOFF-0.5*(A1*FLOAT(N+1)+A2*FLOAT(M+1))
	B1=SAZ*SAL*XSC
	B2=CAZ*SAL*YSC
	B3=ZFAC*CAL
	B4=B3*ZOFF+YOFF-0.5*(B1*FLOAT(N+1)+B2*FLOAT(M+1))
	IAZ=1
	IF (A1.LE.0.0) IAZ=IAZ+1
	IF (A2.LE.0.0) IAZ=IAZ+2
	GOTO (10,20,10,20), IAZ
10	IFIRST=1
	ISTEP=1
	ILAST=M
	GOTO 30
20	IFIRST=M
	ISTEP=-1
	ILAST=1
30	GOTO (50,50,40,40),IAZ
40	JFIRST=1
	JSTEP=1
	JLAST=N
	GOTO 60
50	JFIRST=N
	JSTEP=-1
	JLAST=1
60	GOTO (64,62,62,64),IAZ
62	LLI=1
	GOTO 66
64	LLI=-1
66	IC=ID
	IBEG=IFIRST+ISTEP
70	LNTH=MIN(2*IABS(IBEG-IFIRST)+1,LMAX)
	IF (LLI.EQ.-1) GOTO 72
	LL=0
	GOTO 74
72	LL=LNTH+1
74	I=IBEG
	J=JFIRST
	XX=FLOAT(J)
	YY=FLOAT(I)
	LL=LL+LLI
	D(IX(LL))=A1*XX+A2*YY+A3
	D(IY(LL))=B1*XX+B2*YY+B3*(A(I,J)+ZOFF)+B4
80	I=I-ISTEP
	YY=FLOAT(I)
	LL=LL+LLI
	D(IX(LL))=A1*XX+A2*YY+A3
	D(IY(LL))=B1*XX+B2*YY+B3*(A(I,J)+ZOFF)+B4
	IF (J.EQ.JLAST) GOTO 85
	J=J+JSTEP
	XX=FLOAT(J)
	LL=LL+LLI
	D(IX(LL))=A1*XX+A2*YY+A3
	D(IY(LL))=B1*XX+B2*YY+B3*(A(I,J)+ZOFF)+B4
	IF (I.NE.IFIRST) GOTO 80
85	CALL NXTVU(IC,D,LNTH,W2,NN,IER)
	IF (IER.NE.0) RETURN
	IC=2*ID
	IF (IBEG.EQ.ILAST) GOTO 90
	IBEG=IBEG+ISTEP
	GOTO 70
90	JBEG=JFIRST
100	LNTH=MIN(2*IABS(JBEG-JLAST)+1,LMAX)
	IF (LLI.EQ.-1) GOTO 102
	LL=0
	GOTO 104
102	LL=LNTH+1
104	I=ILAST
	J=JBEG
	XX=FLOAT(J)
	YY=FLOAT(I)
	LL=LL+LLI
	D(IX(LL))=A1*XX+A2*YY+A3
	D(IY(LL))=B1*XX+B2*YY+B3*(A(I,J)+ZOFF)+B4
110	J=J+JSTEP
	XX=FLOAT(J)
	LL=LL+LLI
	D(IX(LL))=A1*XX+A2*YY+A3
	D(IY(LL))=B1*XX+B2*YY+B3*(A(I,J)+ZOFF)+B4
	IF (I.EQ.IFIRST) GOTO 120
	I=I-ISTEP
	YY=FLOAT(I)
	LL=LL+LLI
	D(IX(LL))=A1*XX+A2*YY+A3
	D(IY(LL))=B1*XX+B2*YY+B3*(A(I,J)+ZOFF)+B4
	IF (J.NE.JLAST) GOTO 110
120	CALL NXTVU(ID*2,D,LNTH,W2,NN,IER)
	IF (IER.NE.0) RETURN
	JBEG=JBEG+JSTEP
	IF (JBEG.EQ.JLAST) RETURN
	GOTO 100
500	IER=2
	RETURN
	END
C
C
	FUNCTION ANXTVU(X0,X1,Y0,Y1,X)
C
C	COMPUTE THE EXTENSION OF A LINE FOR NXTVU
C	MODIFIED 8/11/87 BY D.LONG, JPL
C
	IF (X0.EQ.X1) GOTO 10
	ANXTVU=(X-X0)*(Y1-Y0)/(X1-X0)+Y0
	RETURN
10	IF (Y1.GT.Y0) GOTO 20
	ANXTVU=Y0
	RETURN
20	ANXTVU=Y1
	RETURN
	END
C
C
	SUBROUTINE NXT0VU(X,Y,W,KK,LL,IER)
C
C	ROUTINE TO STORE NEW DATA IN SURFACE HISTORY ARRAY IN NXTVU
C	MODIFIED 8/11/87 BY D.LONG, JPL
C
	DIMENSION W(1)
	DATA EPS/0.001/		! ZERO TEST VALUE
	IF (KK.EQ.0) GOTO 10
	IF (2*KK.GE.LL-1) GOTO 20
	IF (ABS(W((KK-1)*2+1)-X)+ABS(W(2*KK)-Y).LT.EPS) RETURN
10	KK=KK+1
	W((KK-1)*2+1)=X
	W(2*KK)=Y
	RETURN
20	IER=1
	RETURN
	END
C
C
	SUBROUTINE HLT3D(A,MD,ND,M,N,W1,L,W2,NN,ALT,AZ,XLEN,XOFF,
     $		YLEN,YOFF,ZFAC,ZOFF,IER)
C
C	ROUTINE TO PLOT COLUMNS IN 3D WITH HIDDEN LINE REMOVAL USING NXTVU
C	WRITTEN 8/17/87 BY D.LONG, JPL
C
C	A	(R): MD BY ND ARRAY CONTAINING Z VALUES TO BE PLOTTED
C	MD,ND	(I): DIMENSIONS OF A
C	M,N	(I): SIZE OF DATA IN A ARRAY  M,N > 1
C	W1	(R): WORKING ARRAY OF DIMENSIONED W1(L)
C		      (NOT USED--INCLUDED FOR COMPATIBILITY WITH PLT3D)
C	L	(I): DIMENSION OF W1 ARRAY
C	W2	(R): WORKING ARRAY OF DIMENSIONED W2(NN)
C	NN	(I): DIMENSION OF W2 ARRAY
C	ALT,AZ	(R): VIEWING ALTITUDE, AZIMUTH IN DEGREES
C	XLEN	(R): LENGTH OF UNPROJECTED X AXIS IN INCHES
C			> 0 SURFACE ABOVE PLANE SHOWN
C			< 0 SURFACE BELOW PLANE SHOWN
C	XOFF	(R): ORIGIN OF PLOT
C	YLEN	(R): LENGTH OF UNPROJECTED Y AXIS IN INCHES
C	YOFF	(R): ORIGIN OF PLOT
C	ZFAC	(R): Z COORDINATE SCALE FACTOR (Zplotted=ZFAC*A(I,J)+ZOFF)
C	ZOFF	(R): Z AXIS OFFSET VALUE
C	IER	(I): RETURNED ERROR CODE
C			0: NO ERROR
C			1: WORKING STORAGE W2 INTERNAL SPACE EXCEEDED
C			2: L LESS THAN 2*MIN(M,N)
C
	DIMENSION A(MD,ND),W2(NN)
	DIMENSION D(12)
C
	COMMON /PLT3B/A1,A2,A3,B1,B2,B3,B4
	DATA DTR/0.0174532925/
C
	IX(I)=(I-1)*2+1
	IY(I)=I*2
C
	TAZ=AZ*DTR
	TALT=ALT*DTR
	SAZ=SIN(TAZ)
	CAZ=COS(TAZ)
	SAL=SIN(TALT)
	CAL=COS(TALT)
	XSC=ABS(XLEN)/FLOAT(N-1)
	ID=1
	IF (XLEN.LT.0.0) ID=-1
	YSC=YLEN/FLOAT(M-1)
	A1=CAZ*XSC
	A2=-SAZ*YSC
	A3=XOFF-0.5*(A1*FLOAT(N+1)+A2*FLOAT(M+1))
	B1=SAZ*SAL*XSC
	B2=CAZ*SAL*YSC
	B3=ZFAC*CAL
	B4=B3*ZOFF+YOFF-0.5*(B1*FLOAT(N+1)+B2*FLOAT(M+1))
	IAZ=1
	IF (A1.LE.0.0) IAZ=IAZ+1
	IF (A2.LE.0.0) IAZ=IAZ+2
	GOTO (10,20,10,20), IAZ
10	IFIRST=1
	ISTEP=1
	ILAST=M
	GOTO 30
20	IFIRST=M
	ISTEP=-1
	ILAST=1
30	GOTO (50,50,40,40),IAZ
40	JFIRST=1
	JSTEP=1
	JLAST=N
	GOTO 60
50	JFIRST=N
	JSTEP=-1
	JLAST=1
60	IC=ID
C
C	PLOT ZOFF PLANE LINE
C
	LL=3
	XX=FLOAT(JFIRST)
	YY=FLOAT(ILAST+ISTEP)+0.001
	D(IX(1))=A1*XX+A2*YY+A3
	D(IY(1))=B1*XX+B2*YY+B3*ZOFF+B4
	XX=FLOAT(JFIRST)
	YY=FLOAT(IFIRST)
	D(IX(2))=A1*XX+A2*YY+A3
	D(IY(2))=B1*XX+B2*YY+B3*ZOFF+B4
	XX=FLOAT(JLAST+JSTEP)+0.001
	YY=FLOAT(IFIRST)
	D(IX(3))=A1*XX+A2*YY+A3
	D(IY(3))=B1*XX+B2*YY+B3*ZOFF+B4
	CALL NXTVU(IC,D,LL,W2,NN,IER)
	IF (IER.NE.0) RETURN
	IC=2*ID
C
	IBEG=IFIRST
	JBEG=JFIRST
70	CONTINUE
C
C	BOTTOM FRONT EDGE
C
	I=IBEG+ISTEP
	J=JBEG
490	CONTINUE
	LL=1
	XX=FLOAT(J)
	YY=FLOAT(I)
	D(IX(1))=A1*XX+A2*YY+A3
	D(IY(1))=B1*XX+B2*YY+B3*ZOFF+B4
	I=I-ISTEP
	LL=LL+1
	YY=FLOAT(I)
	D(IX(2))=A1*XX+A2*YY+A3
	D(IY(2))=B1*XX+B2*YY+B3*ZOFF+B4
	IF (J.EQ.JLAST) GOTO 400
	J=J+JSTEP
	LL=LL+1
	XX=FLOAT(J)
	D(IX(3))=A1*XX+A2*YY+A3
	D(IY(3))=B1*XX+B2*YY+B3*ZOFF+B4
400	CALL NXTVU(IC,D,LL,W2,NN,IER)
	IF (IER.NE.0) RETURN
	IF (I.EQ.IFIRST.OR.J-JSTEP.EQ.JLAST) GOTO 480
	GOTO 490
C
C	PLOT FRONT CENTER LINES
C
480	I=IBEG
	J=JBEG
75	CONTINUE
	XX=FLOAT(J)
	YY=FLOAT(I)
	ZV=ZOFF
	IF (I.NE.IFIRST.AND.J.NE.JFIRST) THEN
		ZV=AMIN1(A(I-ISTEP,J-JSTEP),A(I-ISTEP,J))
		ZV=AMIN1(ZV,A(I,J-JSTEP))+ZOFF
	ENDIF
	D(IX(1))=A1*XX+A2*YY+A3
	D(IY(1))=B1*XX+B2*YY+B3*ZV+B4
	D(IX(2))=D(IX(1))
	D(IY(2))=B1*XX+B2*YY+B3*(A(I,J)+ZOFF)+B4
	CALL NXTVU(IC,D,2,W2,NN,IER)
	IF (IER.NE.0) RETURN
	IC=2*ID
	IF (I.EQ.IFIRST.OR.J.EQ.JLAST) GOTO 80
	I=I-ISTEP
	J=J+JSTEP
	GOTO 75
C
C	FRONT SIDES
C
80	I=IBEG+ISTEP
	J=JBEG
175	CONTINUE
	LL=2
	XX=FLOAT(J)
	YY=FLOAT(I)
	D(IX(1))=A1*XX+A2*YY+A3
	D(IY(1))=B1*XX+B2*YY+B3*ZOFF+B4
	D(IX(2))=D(IX(1))
	D(IY(2))=D(IY(1))+B3*A(I-ISTEP,J)
	IF (I.EQ.IFIRST) GOTO 275
	I=I-ISTEP
	LL=LL+1
	YY=FLOAT(I)
	D(IX(3))=A1*XX+A2*YY+A3
	D(IY(3))=B1*XX+B2*YY+B3*(A(I,J)+ZOFF)+B4
	J=J+JSTEP
	LL=LL+2
	XX=FLOAT(J)
	D(IX(4))=A1*XX+A2*YY+A3
	D(IY(4))=B1*XX+B2*YY+B3*(A(I,J-JSTEP)+ZOFF)+B4
	D(IX(5))=D(IX(4))
	D(IY(5))=D(IY(4))-B3*A(I,J-JSTEP)
275	CALL NXTVU(IC,D,LL,W2,NN,IER)
	IF (IER.NE.0) RETURN
	IF (I.EQ.IFIRST.OR.J-JSTEP.EQ.JLAST) GOTO 260
	GOTO 175
C
C	BACK TOP
C
260	I=IBEG+ISTEP
	J=JBEG
290	CONTINUE
	LL=1
	XX=FLOAT(J)
	YY=FLOAT(I)
	D(IX(1))=A1*XX+A2*YY+A3
	D(IY(1))=B1*XX+B2*YY+B3*(A(I-ISTEP,J)+ZOFF)+B4
	J=J+JSTEP
	LL=LL+1
	XX=FLOAT(J)
	D(IX(2))=A1*XX+A2*YY+A3
	D(IY(2))=B1*XX+B2*YY+B3*(A(I-ISTEP,J-JSTEP)+ZOFF)+B4
	IF (I.EQ.IFIRST) GOTO 300
	I=I-ISTEP
	LL=LL+1
	YY=FLOAT(I)
	D(IX(3))=A1*XX+A2*YY+A3
	D(IY(3))=B1*XX+B2*YY+B3*(A(I,J-JSTEP)+ZOFF)+B4
300	CALL NXTVU(IC,D,LL,W2,NN,IER)
	IF (IER.NE.0) RETURN
	IF (I.EQ.IFIRST.OR.J-JSTEP.EQ.JLAST) GOTO 560
	GOTO 290
C
C	BACK BOTTOM
C
560	I=IBEG+ISTEP
	J=JBEG
590	CONTINUE
	IF (A(I-ISTEP,J).LT.ZOFF) THEN
C
C	FOR COLUMNS BELOW OFFSET PLANE, SHOW BACK CORNER
C
		LL=2
		XX=FLOAT(J+JSTEP)
		YY=FLOAT(I)
		D(IX(1))=A1*XX+A2*YY+A3
		D(IY(1))=B1*XX+B2*YY+B3*(A(I-ISTEP,J)+ZOFF)+B4
		D(IX(2))=A1*XX+A2*YY+A3
		D(IY(2))=B1*XX+B2*YY+B3*ZOFF+B4
		CALL NXTVU(IC,D,LL,W2,NN,IER)
		IF (IER.NE.0) RETURN
	ENDIF
	LL=1
	XX=FLOAT(J)
	YY=FLOAT(I)
	D(IX(1))=A1*XX+A2*YY+A3
	D(IY(1))=B1*XX+B2*YY+B3*ZOFF+B4
	J=J+JSTEP
	LL=LL+1
	XX=FLOAT(J)
	D(IX(2))=A1*XX+A2*YY+A3
	D(IY(2))=B1*XX+B2*YY+B3*ZOFF+B4
	IF (I.EQ.IFIRST) GOTO 550
	I=I-ISTEP
	LL=LL+1
	YY=FLOAT(I)
	D(IX(3))=A1*XX+A2*YY+A3
	D(IY(3))=B1*XX+B2*YY+B3*ZOFF+B4
550	CALL NXTVU(IC,D,LL,W2,NN,IER)
	IF (IER.NE.0) RETURN
	IF (I.EQ.IFIRST.OR.J-JSTEP.EQ.JLAST) GOTO 380
	GOTO 590
C
C	FINISHED ONE LAYER
C
380	CONTINUE
	IF (IBEG.NE.ILAST) THEN
		IBEG=IBEG+ISTEP
		GOTO 70
	ENDIF
	JBEG=JBEG+JSTEP
	IF (JBEG.EQ.JLAST) RETURN	
	GOTO 70
500	IER=2
	RETURN
	END
C
C
	SUBROUTINE NXTVU(IC,D,N,W,NN,IER)
C
C	ROUTINE TO PLOT VISIBLE LINES LINES FOR USE WITH PLT3D
C	LINES BEHIND A PLANAR SURFACE STORED IN W ARE NOT PLOTTED
C	MODIFIED 8/11/87 BY D.LONG, JPL
C
C	IC	(I)	INTIALIZE/DIRECTION CODE
C			 > 0 PLOT SURFACE ABOVE W
C			 < 0 PLOT SURFACE BELOW W
C			 = +1 FIRST CALL FOR ABOVE SURFACE
C			 = +2 FOR SUBSEQUENT CALLS ABOVE SURFACE
C			 = -1 FIRST CALL FOR BELOW SURFACE
C			 NOTE THAT VISIBLE EDGE IS NOT PLOTTED
C			 FOR BELOW PLOT OPTION
C			 = -2 FOR SUBSEQUENT CALLS BELOW SURFACE
C	D	(R)	INPUT POINT PAIRS
C			 D(1) = X(1)
C			 D(2) = Y(1)
C			  ...
C	N	(I)	NUMBER OF POINT PAIRS IN D ARRAY
C	W	(R)	ARRAY STORING TOP (BOTOM) HIDDEN PLANE
C			 (SHOULD NOT BE MODIFIED BETWEEN CALLS)
C	NN	(I)	DIMENSION OF W
C	IER	(I)	RETURNED ERROR FLAG
C			 = 0, OK
C		 	 = 1 OUT OF SPACE ERROR IN W
C
	DIMENSION D(1),W(1)
C
	IX(I)=(I-1)*2+1
	IY(I)=2*I
C
C	DETERMINE DIRECTION AND INITIALIZE
C
	DI=1.0
	IF (IC.LT.0) DI=-1.0
	IF (IABS(IC).NE.1.AND.IC.NE.0) GOTO 20
C
C	CHECK OVERFLOW BUFFER
C
	IF (N.GT.NN/2) GOTO 500
	LL=NN/2-N+1
	I=LL
	W(IX(I))=D(IX(1))
	W(IY(I))=DI*D(IY(1))
	IF (IC.GT.0) CALL PLOT(D(IX(1)),D(IY(1)),3)
	DO 10 J=2,N
		I=I+1
		W(IX(I))=D(IX(J))
		W(IY(I))=DI*D(IY(J))
		IF (IC.GT.0) CALL PLOT(D(IX(J)),D(IY(J)),2)
10	CONTINUE
	IER=0
	RETURN
20	IF (IER.NE.0) RETURN
	II=1
	JJ=LL
	KK=0
	YA0=DI*D(IY(1))
	YB0=W(IY(LL))
	IF (D(IX(1))-W(IX(LL))) 30,30,70
30	CALL PLOT(D(IX(1)),DI*YA0,3)
40	CALL NXT0VU(D(IX(II)),DI*D(IY(II)),W,KK,LL,IER)
	IF (II.EQ.N) GOTO 360
	II=II+1
	YA0=DI*D(IY(II))
	IF (D(IX(II)).GT.W(IX(LL))) GOTO 50
	CALL PLOT(D(IX(II)),DI*YA0,2)
	GOTO 40
50	II=II-1
	XL=D(IX(II))
	YL=DI*D(IY(II))
	YA0=ANXTVU(D(IX(II)),D(IX(II+1)),DI*D(IY(II)),
     1    DI*D(IY(II+1)),W(IX(LL)))
	X0=W(IX(LL))
	IF (YA0.GT.YB0) GOTO 90
	CALL PLOT(X0,DI*YA0,2)
	CALL NXT0VU(X0,YA0,W,KK,LL,IER)
	CALL NXT0VU(X0,YB0,W,KK,LL,IER)
	GOTO 100
70	CALL NXT0VU(W(IX(JJ)),W(IY(JJ)),W,KK,LL,IER)
	IF (JJ.EQ.NN/2) GOTO 380
	JJ=JJ+1
	YB0=W(IY(JJ))
	IF (D(IX(1))-W(IX(JJ))) 80,70,70
80	JJ=JJ-1
	YB0=ANXTVU(W(IX(JJ)),W(IX(JJ+1)),W(IY(JJ)),
     1    W(IY(JJ+1)),D(IX(1)))
	X0=D(IX(1))
	IF (YA0.LE.YB0) GOTO 100
	CALL NXT0VU(X0,YB0,W,KK,LL,IER)
	CALL NXT0VU(X0,YA0,W,KK,LL,IER)
	XL=X0
	YL=YA0
90	IOV0=1
	GOTO 120
100	IOV0=0
120	IF (II.EQ.N) GOTO 300
	IF (JJ.EQ.NN/2) GOTO 310
	IF (D(IX(II+1)).GT.W(IX(JJ+1))) GOTO 130
	ISW=+1
	II=II+1
	X1=D(IX(II))
	YA1=DI*D(IY(II))
	YB1=ANXTVU(W(IX(JJ)),W(IX(JJ+1)),W(IY(JJ)),
     1     W(IY(JJ+1)),X1)
	GOTO 140
130	IF (W(IX(JJ+1)).GT.D(IX(N))) GOTO 340
	ISW=-1
	JJ=JJ+1
	X1=W(IX(JJ))
	YA1=ANXTVU(D(IX(II)),D(IX(II+1)),DI*D(IY(II)),
     1     DI*D(IY(II+1)),X1)
	YB1=W(IY(JJ))
140	IF (YA1.LE.YB1) GOTO 160
	IOV1=1
	IF (IOV0.EQ.0) GOTO 170
150	IF (ISW.EQ.-1) GOTO 200
	CALL NXT0VU(X1,YA1,W,KK,LL,IER)
	CALL PLOT(XL,DI*YL,3)
	CALL PLOT(X1,DI*YA1,2)
	XL=X1
	YL=YA1
	GOTO 200
160	IOV1=0
	IF (IOV0.EQ.0) GOTO 190
170	CONTINUE
	IF (YA1-YB1+YB0-YA0.EQ.0.0) THEN
		FRAC=1.0
	ELSE
		FRAC=(YB0-YA0)/(YA1-YB1+YB0-YA0)
	ENDIF
	XI=(X1-X0)*FRAC+X0
	YI=(YA1-YA0)*FRAC+YA0
	CALL NXT0VU(XI,YI,W,KK,LL,IER)
	IF (IOV0.EQ.0) GOTO 180
	CALL PLOT(XL,DI*YL,3)
	CALL PLOT(XI,DI*YI,2)
	XL=XI
	YL=YI
	GOTO 190
180	XL=XI
	YL=YI
	GOTO 150
190	IF (ISW.EQ.+1) GOTO 200
	CALL NXT0VU(W(IX(JJ)),W(IY(JJ)),W,KK,LL,IER)
200	IF (IER.NE.0) RETURN
	X0=X1
	YA0=YA1
	YB0=YB1
	IOV0=IOV1
	GOTO 120
310	X1=W(IX(NN/2))
	YA1=ANXTVU(D(IX(II)),D(IX(II+1)),DI*D(IY(II)),
     1     DI*D(IY(II+1)),X1)
	YB1=W(IY(NN/2))
	IF (YA1.GT.YB1) GOTO 320
	CALL NXT0VU(X1,YB1,W,KK,LL,IER)
	CALL NXT0VU(X1,YA1,W,KK,LL,IER)
	CALL PLOT(X1,DI*YA1,3)
	GOTO 330
380	II=1
320	CALL PLOT(D(IX(II)),D(IY(II)),3)
330	IF (II.EQ.N) GOTO 400
	II=II+1
	CALL NXT0VU(D(IX(II)),DI*D(IY(II)),W,KK,LL,IER)
	CALL PLOT(D(IX(II)),D(IY(II)),2)
	GOTO 330
300	IF (JJ.EQ.NN/2) GOTO 400
340	X1=D(IX(N))
	YA1=DI*D(IY(N))
	YB1=ANXTVU(W(IX(JJ)),W(IX(JJ+1)),W(IY(JJ)),
     1     W(IY(JJ+1)),X1)
	IF (YA1.LE.YB1) GOTO 350
	CALL NXT0VU(X1,YA1,W,KK,LL,IER)
	CALL NXT0VU(X1,YB1,W,KK,LL,IER)
	CALL PLOT(X1,DI*YA1,2)
350	IF (JJ.EQ.NN/2) GOTO 400
	JJ=JJ+1
	CALL NXT0VU(W(IX(JJ)),W(IY(JJ)),W,KK,LL,IER)
	GOTO 350
360	JJ=0
	GOTO 350
400	LL=NN/2-KK+1
	I=LL
	DO 410 J=1,KK
		W(IX(I))=W(IX(J))
		W(IY(I))=W(IY(J))
		I=I+1
410	CONTINUE
	RETURN
500	IER=1
	RETURN
	END
